A random wind field fast simulation method and system based on a double-layer parallel architecture

Abstract

The application discloses a random wind field fast simulation method and system based on a double-layer parallel architecture, and the method is based on a Shinozuka harmonic synthesis framework, adopts a space+frequency double-layer parallel architecture, and removes nested loops through a tensor broadcast mechanism in a space dimension, so that a coherent weighted distance matrix of all space point pairs is obtained through single parallel calculation; the single-frequency amplitude calculation is encapsulated as a pure function in a frequency dimension, a multi-frequency task batch parallel is realized through a vectorization mapping operator, and a self-adaptive memory management mechanism is matched to dynamically adjust the batch size. The system is equipped with parameter input, parallel calculation, memory management, time domain generation and other modules, and supports CPU / GPU heterogeneous acceleration. The application significantly improves the random wind field simulation efficiency on the premise of maintaining complete physical precision, and is suitable for wind resistance design of large-span bridges, super high-rise buildings, offshore wind turbine units and the like.

Description

Technical Field

[0001] This invention belongs to the field of electrical digital data processing technology, and relates to a parallel acceleration method for multi-dimensional batch numerical computation and a heterogeneous computing resource scheduling system, especially to a fast simulation method and system for stochastic wind fields based on a two-layer parallel architecture. Background Technology

[0002] With the popularization of CPU / GPU heterogeneous computing hardware, the demand for multi-dimensional batch tensor operations in scenarios such as stochastic process simulation, signal processing, and scientific computing is growing rapidly. These scenarios generally involve matrix operations coupled in two dimensions of space and frequency, as well as batch numerical superposition calculations, which have extremely high requirements for hardware parallel resource utilization and memory scheduling stability.

[0003] Current general parallel solutions for tensor operations in the field of electronic digital data processing generally suffer from three common technical defects: First, they often use nested loops to traverse computational units of different dimensions, resulting in a utilization rate of less than 30% for CPU SIMD instruction sets and GPU CUDA cores, leading to serious waste of computing power; second, they only support single-dimensional parallel scheduling, and when the number of computational units increases significantly, the computation time increases linearly, making it unsuitable for large-scale scenarios; and third, they adopt fixed-batch task scheduling strategies, failing to dynamically adapt to the remaining memory on the device, resulting in a memory overflow probability exceeding 15% during computation and poor operational stability.

[0004] The Shinozuka harmonic synthesis framework is the most widely used technique in the field of stochastic process simulation. Its computational process involves coupled batch tensor operations in both spatial and frequency dimensions, and stochastic wind field simulation is a typical application scenario for this framework. As civil engineering structures become larger and more flexible, structures such as long-span bridges, super high-rise buildings, and offshore wind turbines are becoming increasingly sensitive to wind loads. Accurate simulation of spatially correlated stochastic wind fields is fundamental to the wind-resistant design and safety assessment of such structures. The simulation results are key inputs for aeroelastic analysis, buffeting response prediction, fatigue assessment, and applicability evaluation. With the increasing scale and complexity of structures, the need for large-scale spatial wind field simulations involving hundreds or even thousands of monitoring points is becoming increasingly urgent.

[0005] When simulating random wind fields using the classic Shinozuka harmonic synthesis method, the following steps need to be performed: (1) Calculate the autospectral density of each spatial point based on the wind spectrum model (such as the Kaimal spectrum); (2) Construct the cross-spectral density matrix by combining coherence function models (such as Davenport exponential decay); (3) Perform Cholesky decomposition on the cross spectral density matrix; (4) Combined with a random phase, it generates a complex amplitude; (5) Synthesize the time-domain wind speed sequence by inverse Fourier transform.

[0006] This method requires constructing an n×n dimensional cross-spectral density matrix at each frequency point and iteratively processing all frequencies, resulting in a computational complexity of O(n^2). 2 When the number of spatial points n or frequency bands N is large, there are serious computational and memory bottlenecks. Traditional methods take hours or even days, which is difficult to meet the needs of large-scale simulation in modern engineering. Existing wind field simulation acceleration methods mostly use interpolation or low-rank matrix decomposition techniques. By reducing frequency segments or using low-rank matrices to approximate the original samples, computational efficiency can be improved, but simulation accuracy is sacrificed to some extent.

[0007] Therefore, there is an urgent need for a stochastic wind field simulation technology that can maintain complete physical accuracy, accurately reproduce the target spectral characteristics and spatial correlation structure, and significantly improve computational efficiency and make full use of modern parallel computing hardware. Summary of the Invention

[0008] To address the aforementioned technical problems, the first objective of this invention is to provide a fast simulation method for stochastic wind fields based on a two-layer parallel architecture, and the second objective is to provide a fast simulation system for stochastic wind fields based on a two-layer parallel architecture.

[0009] To achieve the aforementioned primary objective, this invention provides a technical solution for a fast simulation method of stochastic wind fields based on a dual-layer parallel architecture. The dual-layer parallelism includes two layers: spatial dimension parallelism and frequency dimension parallelism. The method is implemented based on a harmonic superposition-type stochastic process simulation framework, and includes: To address the need for calculating the spatial correlation of wind fields at multiple spatial points, a tensor broadcasting mechanism is adopted to achieve element-level parallel computation: the spatial coordinate vectors are reshaped into n×1 and 1×n dimensional tensors respectively, and the numerical computation framework automatically expands and matches the dimensions. A single parallel computation yields the spatial coherence weighted distance matrix of all spatial point pairs, where n is the total number of spatial points. The calculation of wind field amplitude at a single frequency is encapsulated into a pure function with no side effects that only takes the frequency value and random phase vector as input and outputs the amplitude vector. The pure function is then frequency vectorized by a vectorization mapping operator, and the amplitude calculation tasks of multiple frequencies are executed in batches to achieve parallelism in the frequency dimension.

[0010] Furthermore, it includes the following steps: S1: Obtain simulation parameters, including: number of spatial points n, number of frequency bands N, number of time points M, and spatial coordinates. Average wind speed vector Angular frequency vector Random phase matrix Angular frequency increment Δω, spatial coherence attenuation coefficient Surface friction wind speed u* Where the angular frequency ω is related to the physical frequency n f Satisfying ω=2πn f ; S2: Parallel computation of the spatial coherence weighted distance matrix in spatial dimensions: The spatial coordinate vectors x, y, z are reshaped into n×1 and 1×n tensors, respectively, denoted as... and Using the tensor broadcasting mechanism, the matrix is ​​automatically expanded to an n×n dimension through element-level operations, and the spatial coherence weighted distance matrix is ​​obtained through parallel computation. The calculation formula is as follows: Formula (1) In formula (1), : an n×n dimensional spatial coherent weighted distance matrix, where D i,j This represents the distance between the i-th spatial point and the j-th spatial point; Spatial coherence attenuation coefficients in the x, y, and z directions are used to describe the degree of attenuation of wind speed correlation in different directions; : The x, y, z coordinates of the i-th spatial point, with a tensor dimension of n×1; : The x, y, z coordinates of the j-th spatial point, with a tensor dimension of 1×n; The tensor broadcasting mechanism refers to the mechanism by which the numerical computing framework automatically expands the dimensions of tensors of different dimensions to match the shape of the tensors in order to support element-level operations, without the need to manually write loop logic. S3: Define a single-frequency amplitude vector B (n) The spatially parallel computation function single_freq_amplitude includes the following sub-steps: S31: Transform the average wind speed vector U z Remodeled into n×1 and 1×n dimensional tensors, denoted as U zi and U zj ; S32: Calculate the coherence function matrix: Formula (2) In formula (2), : an n×n dimensional Davenport coherence function matrix, where Γ i,j This represents the coherence coefficient between the i-th spatial point and the j-th spatial point; ω: The currently calculated angular frequency value; D: Spatial coherence weighted distance matrix, which is the same as D in formula (1). (n×n) same; U zi: The average wind speed at the i-th spatial point, with a tensor dimension of n×1; U zj : The average wind speed at the j-th spatial point, with a tensor dimension of 1×n; S33: Calculate the autospectral density of each spatial point based on the Davenport downwind wind spectrum model, and reshape it into n×1 and 1×n dimensional tensors, denoted as S. i and S j ; S34: Calculate the cross-spectral density matrix: Formula (3) In formula (3), : an n×n dimensional cross-spectral density matrix, where S i,j This represents the cross-spectral density between the i-th spatial point and the j-th spatial point; : Davenport coherence function matrix, and Γ in formula (2) (n×n) same; ⊙: Element-wise multiplication, multiplying corresponding elements; S i : The autospectral density of the i-th spatial point, with tensor dimension n×1; : The autospectral density of the j-th spatial point, where the superscript T is the matrix transpose operator and the tensor dimension is 1×n; S35: Perform Cholesky decomposition on the cross-spectral density matrix to obtain the lower triangular matrix H. (n×n) So that S=HH T ; S36: Calculate the amplitude vector at the current frequency: Formula (4) In formula (4), : n-dimensional amplitude vector, where B i This represents the amplitude of the i-th spatial point at the current frequency; An n×n lower triangular matrix, obtained by Cholesky decomposition of the cross-spectral density matrix, satisfying S=HH T ; i: Imaginary unit, i 2 =-1; : an n-dimensional random phase vector, where This represents the random phase of the i-th spatial point at the current frequency, with a value range of [0, 2π]. S4: Adaptive memory management: Dynamically determines the frequency and batch size, including: First, estimate the memory requirements for frequency calculation in a single batch, and then calculate the optimal batch size N based on the available device memory. batch ; Then divide the N frequencies into One batch, This is for rounding up; S5: Parallel computation of the complete amplitude matrix in the frequency dimension B (N×n) The single_freq_amplitude function defined in step S3 is vectorized using the vectorized mapping operator vmap, with a batch size of N for each calculation. batch The frequency is calculated, and after K batch calculations are performed, the complete amplitude matrix is ​​obtained by splicing them together; S6: Generate time-domain wind speed sequence: Perform transpose, zero-filling, inverse fast Fourier transform, and phase modulation on the frequency domain amplitude matrix to generate the time-domain wind speed sequence f. (n×M) The number of time points M satisfies M≥2N to avoid frequency domain aliasing.

[0011] Furthermore, in step S35, a positive semidefinite correction is performed on the cross-spectral density matrix before Cholesky decomposition to avoid task interruption due to decomposition failure. The positive semidefinite correction rule is as follows: When an accuracy better than 0.1% is required, the disturbance coefficient σ is taken as 10. -8 For conventional engineering calculations, σ = 10. -7 When the allowable accuracy error is ≥1%, take σ=10. -6 The correction method involves adding a small perturbation term σI to the cross-spectral density matrix, where I is the identity matrix; if σ=10 is used... -6 If the Cholesky decomposition still fails, the input parameter validation process is automatically triggered, returning the location of the abnormal parameters and modification prompts.

[0012] Furthermore, step S4 includes the following sub-steps: S41: Single Batch Memory Requirements Calculation: Formula (5) In formula (5), : Memory required for a single batch frequency calculation, in GB; n: number of space points, 3n 2 The memory usage corresponds to three types of n×n-dimensional tensors: spatial coherence weighted distance matrix, coherence function matrix, and cross-spectral density matrix. 4n corresponds to four types of n-dimensional tensors: average wind speed vector, self-spectral density vector, random phase vector, and amplitude vector. The number of bytes for a floating-point data type, such as 4 bytes for a 32-bit single-precision floating-point number and 8 bytes for a 64-bit double-precision floating-point number; : Frequency of processing per batch; S42: Set safety factor ,judge Is it less than or equal to the available memory M? avail ,in For safety factors, the value ranges from 1.2 to 1.5; S43: If the limit is exceeded, calculate the optimal batch size using the following formula: Formula (6) In formula (6), Optimal batch size, i.e., the frequency of processing each time; : Round down; Available memory size, in GB; Safety factor, ranging from 1.2 to 1.5, is used to reserve a certain amount of memory space to avoid memory overflow; S44: Divide the N frequencies into One batch, for use in subsequent calculations.

[0013] Furthermore, the calculation method for the complete amplitude matrix in step S5 includes first calculating... The amplitude matrix is ​​transposed and then concatenated to form the complete amplitude matrix. : The The formula for calculating the amplitude matrix is: Formula (7) In formula (7), : The amplitude matrix is ​​dimensional, where B k,i This represents the amplitude of the i-th spatial point at the k-th frequency; Vectorized mapping operator, used to apply a single-frequency amplitude calculation function in batches to multiple frequencies; : Single-frequency amplitude calculation function, the input is frequency value and random phase vector, and the output is amplitude vector; :N batch A 3D angular frequency vector containing all frequency values ​​processed in the current batch; : A dimensional random phase matrix, where Φ k,iThis represents the random phase of the i-th spatial point at the k-th frequency.

[0014] Will 3D amplitude matrix The transpose operation is ; The formula for splicing the complete amplitude matrix is: Formula (8) In formula (8), : an n×N dimensional complete amplitude matrix, containing the amplitude of all spatial points at all frequencies; The amplitude matrix of the kth batch, with dimensions n×N batch ; K: Total number of batches ,in This is for rounding up.

[0015] Furthermore, the zero-filling, inverse fast Fourier transform, and phase modulation in step S6 include: S61: Zero-filling: Concatenate an n×(MN) dimensional zero matrix to obtain the zero-filled amplitude matrix: Formula (9) In formula (9), The amplitude matrix after zero-filling has a dimension of n×M; : The original amplitude matrix, with dimensions n×N; A zero matrix with n rows (MN) and columns, used to extend the number of columns of the amplitude matrix from N to M; M: Number of time points, i.e., the length of the generated time-domain wind speed sequence; S62: Inverse Fast Fourier Transform, the calculation formula is: Formula (10) In formula (10), The complex matrix after inverse Fourier transform has dimensions n×M; IFFT(·): Inverse Fast Fourier Transform operator, used to convert frequency domain signals into time domain signals; The amplitude matrix after zero-filling is the same as B in formula (9). padded same; S63: Phase modulation, calculated using the following formula: Formula (11) In formula (11), : n×M dimensional time-domain wind speed sequence, where f i,t This represents the wind speed value at the i-th spatial point at the t-th time point; Angular frequency increment, Δω=ω up / N, where ω up The upper limit frequency; Re{·}: The operation of taking the real part of a complex number; The complex matrix after the inverse Fourier transform, and G in formula (10) (n×M) same; Time point index vector ; i: Imaginary unit, i 2 =-1.

[0016] Further, in step S33, the formula for calculating the Davenport downwind wind spectrum is: Formula (12) In formula (12), : Power spectral density of wind speed in the downwind direction; u * : Surface friction wind speed, in m / s; n f Physical frequency, measured in Hz; f: dimensionless frequency , where U is the average wind speed at height z in m / s, and z is the height above the ground in m.

[0017] Furthermore, the implementation of the vectorized mapping operator vmap in step S5 includes: (a) JAX backend: Calls the jax.vmap function and combines it with jax.jit for just-in-time compilation optimization; (b) PyTorch backend: Calls the torch.vmap function and uses CUDA for GPU acceleration.

[0018] Furthermore, in step S5, the amplitude matrix calculation for each batch of frequencies is performed in parallel, with no data dependency between batches. Instead, a serial processing method is used. After each batch is processed, the memory occupied by intermediate variables, including the cross-spectral density matrix and the lower triangular matrix of Cholesky decomposition, is released, effectively reducing the overall memory usage.

[0019] Furthermore, the inverse fast Fourier transform described in step S6 adopts a fully parallel computing method, performing inverse fast Fourier transform operations simultaneously on the frequency domain sequences of n spatial points.

[0020] Furthermore, if a JAX computing backend is used, step S0 is also included: before executing step S1, just-in-time (JIT) optimization is performed on the single-frequency amplitude calculation function and the vectorized mapping operator vmap to compile the function into efficient machine code, thereby improving the overall computing performance.

[0021] To achieve the second objective mentioned above, the present invention provides a technical solution for a fast simulation system of stochastic wind fields based on a two-layer parallel architecture, comprising: (a) Parameter input module: used to acquire and store parameters required for simulation, including spatial coordinates, wind speed, frequency, phase, and coherence coefficient; (b) Spatial distance matrix calculation module: used to reshape spatial coordinates into tensors, achieve dimension matching through tensor broadcasting mechanism, calculate spatial coherent weighted distance matrix in parallel, and output the distance matrix of all spatial point pairs in one go without nested loops; (c) Single-frequency amplitude calculation module: Defines the spatial parallel calculation function for single-frequency amplitude vector, including steps for coherence function matrix calculation, cross-spectral density matrix calculation, Cholesky decomposition and amplitude vector calculation. The entire process adopts tensor broadcasting mechanism to achieve spatial parallelism. (d) Parameter verification and matrix correction module: used to perform positive semidefinite correction on the cross-spectral density matrix, and to trigger the input parameter verification process when the decomposition fails, returning the location of abnormal parameters and modification prompts; (e) Adaptive memory management module: including memory reading unit, memory estimation unit and adaptive batching unit. The memory reading unit is used to obtain the available memory of the current device through the operating system API. The memory estimation unit is used to add a frame correction factor to estimate the memory requirement of a single batch. The adaptive batching unit is used to calculate the optimal batch size and process the batches when the memory requirement exceeds the limit. (f) Frequency parallel amplitude matrix calculation module: Using the vectorized mapping operator vmap, the amplitude matrix of a frequency batch is calculated each time, and finally spliced ​​to obtain the complete amplitude matrix; (g) Time-domain wind speed calculation module: used to perform zero-filling, inverse fast Fourier transform and phase modulation on the frequency domain amplitude matrix to generate a time-domain wind field sequence.

[0022] Furthermore, it also includes a JIT compilation module: used in the JAX backend to perform just-in-time compilation optimization by setting static parameters for the single-frequency amplitude calculation function and the vectorized mapping operator vmap.

[0023] To achieve the above objectives, the present invention provides a technical solution for an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the fast simulation method for stochastic wind fields based on a two-layer parallel architecture.

[0024] To achieve the above objectives, the present invention provides a computer-readable storage medium storing computer program instructions, wherein the computer program, when executed by a processor, implements the fast simulation method for stochastic wind fields based on a two-layer parallel architecture.

[0025] Compared with existing technologies, the beneficial effects achieved by this invention are as follows: 1. Significantly improves computational efficiency: Through dual-layer parallel computing and vectorized mapping technology, the computation time for stochastic wind field simulation is greatly shortened, making it suitable for simulation tasks with a large number of spatial points and frequency bands; 2. Significantly optimized stability of heterogeneous computing resource scheduling: The adaptive memory management mechanism can dynamically read the available memory of the device through the operating system API and automatically adjust the computing batch size, completely solving the memory overflow problem in high computing power scenarios. It is compatible with all levels of hardware environment from embedded devices to high-performance computing clusters and can flexibly adapt to computing tasks of different scales.

[0026] 3. Maintain complete physical accuracy: It adopts the same physical model and mathematical foundation as the traditional Shinozuka method, avoiding the accuracy loss caused by interpolation or low-rank approximation, and accurately reproduces the target wind spectrum characteristics and spatial correlation structure; 4. Flexible adaptation to different problem scales: The adaptive memory management mechanism ensures that the method of this invention can run stably under computational tasks of different scales, and is suitable for large-scale engineering simulations ranging from hundreds to thousands of spatial points; 5. Strong adaptability to multiple hardware / frameworks: It natively supports two mainstream numerical computing frameworks, JAX and PyTorch, and automatically adapts to heterogeneous computing resources such as CPU / GPU. It can achieve cross-platform operation without requiring users to manually modify the code, reducing the adaptation cost of technology implementation. 6. Wide range of engineering applications: Applicable to wind field simulation of various civil engineering structures such as long-span bridges, super high-rise buildings, and offshore wind turbines, providing an efficient and reliable calculation tool for structural wind resistance design and safety assessment. It can reduce the time required for large-scale wind field simulation from several hours to seconds, greatly improving the efficiency of engineering simulation. Attached Figure Description

[0027] Figure 1 This is a flowchart illustrating a fast simulation method for stochastic wind fields based on a two-layer parallel architecture, according to the present invention.

[0028] Figure 2This is a comparison chart of simulated wind speed time series at different spatial locations at different altitudes in Embodiment 1 of the present invention.

[0029] Figure 3 This is a schematic diagram of the distribution of ADF stationarity test p-values ​​for all spatial point simulated wind speed sequences in Embodiment 1 of the present invention.

[0030] Figure 4 In the figure, 'af' is a schematic diagram comparing the simulated wind speed power spectral density at different spatial locations with the Davenport theoretical spectrum in Embodiment 1 of the present invention.

[0031] Figure 5 In the diagram, 'ad' represents a comparison between the simulated cross-correlation function and the theoretical correlation function for different spatial locations in Embodiment 1 of the present invention.

[0032] Figure 6a The graph shows the comparison of the total computation time for different spatial points at a fixed frequency range of N=3000 in the computation efficiency test results of Embodiment 2 of the present invention.

[0033] Figure 6b This is a comparison curve of the total calculation time for different frequency bands in the calculation efficiency test results of Embodiment 2 of the present invention, with a fixed number of spatial points n=100.

[0034] Figure 4 The af contains 6 sub-maps, each corresponding to a specific spatial point. The vertical axis of the sub-map is the frequency-weighted normalized downwind power spectral density, which is obtained by multiplying the output of formula (12) by the physical frequency n. f The horizontal axis is the dimensionless frequency f defined by formula (12), which shows the comparison between the power spectral density (blue curve) of the simulated wind speed at the sub-plot location and the theoretical spectrum (orange curve) under frequency weighted normalization. This is used to verify whether the spectral characteristics of the simulated wind speed are consistent with the theory.

[0035] Figure 5 The ad plot contains four subplots, each showing the comparison between the simulated wind speed cross-correlation function (blue curve) and the theoretical correlation function (orange curve) for a specific spatial point pair at different time lags. This is used to verify whether the spatial correlation structure of the simulated wind speed is consistent with the theory. Detailed Implementation

[0036] To make the above-mentioned objectives, features, and advantages of this application more apparent and understandable, the specific embodiments of this application will be described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the scope of this application. Furthermore, it should be noted that, for ease of description, only the parts relevant to this application are shown in the accompanying drawings, not the entire structure. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without inventive effort are within the scope of protection of this application.

[0037] In this document, the term "embodiment" means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of this application. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a separate or alternative embodiment mutually exclusive with other embodiments. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0038] This invention proposes a fast simulation method for stochastic wind fields based on a two-layer parallel architecture. All steps can be executed by an electronic device containing a processor and memory. It is applicable to numerical simulation scenarios of stochastic processes based on the Shinozuka harmonic synthesis framework. This embodiment uses stochastic wind field simulation as a typical application scenario for illustration. The core method can be reused to all space-frequency two-dimensional coupled batch tensor operation scenarios. Figure 1 As shown, it includes the following steps:

[0039] S1: Acquire and verify simulation parameters. The processor receives simulation parameters through the parameter input interface, including: number of spatial points n, number of frequency bands N, number of time points M, and spatial coordinates. Average wind speed vector Angular frequency vector Random phase matrix Angular frequency increment Δω, spatial coherence attenuation coefficient Surface friction wind speed u * Where the angular frequency ω is related to the physical frequency n f Satisfying ω=2πn f ; The processor automatically performs parameter validity checks: it checks whether the average wind speed is all positive, whether there are duplicate points in the spatial coordinates, and whether the coherence coefficient is within the reasonable range of [0,50]. If any abnormality is found, a corresponding prompt will be displayed.

[0040] S2: Spatial Dimension Parallel Computation of Distance Matrix: The processor invokes the tensor broadcast operation channel of the numerical computation framework to reshape the spatial coordinate vectors x, y, z into n×1 and 1×n dimensional tensors, respectively, denoted as... and ; Element-level parallel computation is performed using the CPU SIMD instruction set or the GPU CUDA core, automatically expanding to an n×n dimensional matrix, and the coherent weighted distance matrix of all spatial point pairs is obtained through parallel computation: Formula (1) In formula (1), : an n×n dimensional spatial coherent weighted distance matrix, where D i,j This represents the distance between the i-th spatial point and the j-th spatial point; Spatial coherence attenuation coefficients in the x, y, and z directions are used to describe the degree of attenuation of wind speed correlation in different directions; : The x, y, z coordinates of the i-th spatial point, with a tensor dimension of n×1; : The x, y, z coordinates of the j-th spatial point, with a tensor dimension of 1×n; Compared to the traditional nested loops (see Algorithm 1 below), this step reduces the computational complexity from O(n^2) to O(n^2). 2 The computation is reduced from ×N) to a single vectorized operation, which can be parallelized on the CPU using SIMD instructions, and thousands of threads can be invoked on the GPU for simultaneous computation.

[0041] S3: Single-frequency amplitude function encapsulation: Define single-frequency amplitude vector B (n) The spatial parallel computation function single_freq_amplitude encapsulates the amplitude calculation at a single frequency into a pure function without side effects. It only takes the current frequency value and the random phase vector of the current frequency as input and outputs the amplitude vector of that frequency. The intermediate variables (coherence matrix, cross-spectral matrix, and lower triangular matrix of Cholesky decomposition) are automatically released by the processor after calculation, without the need for manual management. All operations inside the function use the tensor broadcast mechanism to achieve spatial parallelism. The calculation process includes calculation formulas (2)-(4). The processor automatically performs positive semidefinite correction on the cross-spectral density matrix to avoid Cholesky decomposition failure.

[0042] Includes the following sub-steps: S31: Transform the average wind speed vector U z Remodeled into n×1 and 1×n dimensional tensors, denoted as U zi and U zj Automatic broadcast matching dimensions; S32: Calculation of coherence function matrix: Formula (2) In formula (2), : an n×n dimensional Davenport coherence function matrix, where Γi,j This represents the coherence coefficient between the i-th spatial point and the j-th spatial point; ω: The currently calculated frequency value; D: Spatial coherence weighted distance matrix, which is the same as D in formula (1). (n×n) same; U zi : The average wind speed at the i-th spatial point, with a tensor dimension of n×1; U zj : The average wind speed at the j-th spatial point, with a tensor dimension of 1×n.

[0043] S33: Based on the empirical wind spectrum model, calculate the autospectral density of each spatial point and reshape it into n×1 and 1×n dimensional tensors, denoted as S. i and S j ; The empirical wind spectrum model is the Davenport downwind wind spectrum, and the calculation formula is: Formula (12) In formula (12), : Power spectral density of wind speed in the downwind direction; u * : Surface friction wind speed, in m / s; n f Physical frequency, measured in Hz; f: dimensionless frequency , where U is the average wind speed at height z in m / s, and z is the height above the ground in m.

[0044] S34: Calculate the cross-spectral density matrix: Formula (3) In formula (3), : an n×n dimensional cross-spectral density matrix, where S i,j This represents the cross-spectral density between the i-th spatial point and the j-th spatial point; : Davenport coherence function matrix, and Γ in formula (2) (n×n) same; ⊙: Element-wise multiplication, multiplying corresponding elements; S i : The autospectral density of the i-th spatial point, with tensor dimension n×1; S j : The autospectral density of the j-th spatial point, with a tensor dimension of 1×n.

[0045] S35: Perform Cholesky decomposition on the cross-spectral density matrix to obtain the lower triangular matrix H. (n×n) So that S=HH T ; The cross-spectral density matrix obtained during the calculation is subjected to positive semidefinite correction before Cholesky decomposition to avoid task interruption due to decomposition failure. The positive semidefinite correction rule is as follows: A small perturbation term σI is added to the cross-spectral density matrix, where I is the identity matrix. The value of σ is determined based on the required accuracy: σ = 10 when a calculation accuracy better than 0.1% is required. -8 For conventional engineering calculations, σ = 10. -7 When the allowable accuracy error is ≥1%, take σ=10. -6 ; If σ is 10 -6 If the Cholesky decomposition still fails, the input parameter verification process is automatically triggered, prompting the user to check the rationality of the spatial coordinates and average wind speed parameters.

[0046] S36: Calculate the amplitude vector: Formula (4) In formula (4), : n-dimensional amplitude vector, where B i This represents the amplitude of the i-th spatial point at the current frequency; An n×n lower triangular matrix, obtained by Cholesky decomposition of the cross-spectral density matrix, satisfying S=HH T ; i: Imaginary unit, i 2 =-1; : an n-dimensional random phase vector, where It represents the random phase of the i-th spatial point at the current frequency, with a value range of [0, 2π].

[0047] S4: Adaptive Memory Management: The processor reads the available memory size of the current device in real time through the operating system API, estimates the memory requirements for single-batch frequency calculation according to formula (5), and dynamically adjusts the segment size of frequency batch processing in combination with the preset safety factor to ensure that no memory overflow occurs during frequency parallel calculation. This includes the following steps: S41: Estimate the memory requirements for a single batch frequency calculation: Formula (5) In formula (5), : Memory required for a single batch frequency calculation, in GB; n: number of space points, 3n 2 The memory usage corresponds to three types of n×n-dimensional tensors: spatial coherence weighted distance matrix, coherence function matrix, and cross-spectral density matrix. 4n corresponds to four types of n-dimensional tensors: average wind speed vector, self-spectral density vector, random phase vector, and amplitude vector. The number of bytes for a floating-point data type, such as 4 bytes for a 32-bit single-precision floating-point number and 8 bytes for a 64-bit double-precision floating-point number; : Frequency of processing per batch; S42: Judgment Is it less than or equal to the available memory M? avail ,in For safety margin, it is generally taken as 1.2 to 1.5; S43: If Then adjust downwards. Until the memory constraint is met, the optimal batch size is determined. The calculation formula is: Formula (6) In formula (6), Optimal batch size, i.e., the frequency of processing each time; : Round down; Available memory size, in GB; Safety factor, typically between 1.2 and 1.5, is used to reserve a certain amount of memory space to avoid memory overflow.

[0048] S44: Divide the N frequencies into One batch, for use in subsequent calculations; S5: Parallel computation of the complete amplitude matrix using frequency B (N×n) The processor calls the vectorization mapping operator vmap to perform frequency vectorization mapping on the single_freq_amplitude function defined in step S3. It supports two computing backends: JAX and PyTorch. The JAX backend calls jax.vmap and combines it with JIT just-in-time compilation optimization, while the PyTorch backend calls torch.vmap and uses CUDA for GPU acceleration. The intermediate variable memory is automatically released after each batch of calculations is completed. After all batches of calculations are completed, the complete amplitude matrix is ​​obtained by concatenating along the frequency dimension. The calculation of the complete amplitude matrix in step S5 includes first calculating... The amplitude matrix is ​​transposed and then concatenated to form the complete amplitude matrix. : The The formula for calculating the amplitude matrix is: Formula (7) In formula (7), : The amplitude matrix is ​​dimensional, where B k,i This represents the amplitude of the i-th spatial point at the k-th frequency; Vectorized mapping operator, used to apply a single-frequency amplitude calculation function in batches to multiple frequencies; : Single-frequency amplitude calculation function, the input is frequency value and random phase vector, and the output is amplitude vector; :N batch A 3D angular frequency vector containing all frequency values ​​processed in the current batch; : A dimensional random phase matrix, where Φ k,i This represents the random phase of the i-th spatial point at the k-th frequency; The vectorized mapping operator refers to a general parallel operator that automatically applies single-input, single-output functions to batch inputs. Different computing frameworks have corresponding implementations (jax.vmap in JAX, torch.vmap in PyTorch, etc.).

[0049] Will 3D amplitude matrix The transpose operation is ; The formula for splicing the complete amplitude matrix is: Formula (8) In formula (8), : an n×N dimensional complete amplitude matrix, containing the amplitude of all spatial points at all frequencies; The amplitude matrix of the kth batch, with dimensions n×N batch ; K: Total number of batches ,in This is for rounding up.

[0050] S6: Generating the Time-Domain Wind Speed ​​Sequence: The processor performs zero-filling, inverse fast Fourier transform, and phase modulation operations on the complete amplitude matrix, and simultaneously performs parallel inverse Fourier transforms on the frequency domain sequences of n spatial points, finally outputting a pulsating wind field time-domain wind speed sequence f that conforms to the target wind spectrum characteristics. (n×M) This includes the following steps: Step S6 includes the following sub-steps: S61: Zero-filling: Concatenate an n×(MN) dimensional zero matrix to obtain the zero-filled amplitude matrix: Formula (9) In formula (9), The amplitude matrix after zero-filling has a dimension of n×M; : The original amplitude matrix, with dimensions n×N; : A zero-dimensional matrix is ​​used to extend the number of columns in the amplitude matrix from N to M. M: Number of time points, i.e., the length of the generated time-domain wind speed sequence; S62: Inverse Fast Fourier Transform: Perform the inverse fast fourier transform simultaneously on the frequency domain sequence of all spatial points to obtain a complex time domain matrix; The formula for the inverse fast Fourier transform is: Formula (10) In formula (10), The complex matrix after inverse Fourier transform has dimensions n×M; IFFT(·): Inverse Fast Fourier Transform operator, used to convert frequency domain signals into time domain signals; The amplitude matrix after zero-filling is the same as B in formula (9). padded same; S63: Phase Modulation: Formula (11) In formula (11), : n×M dimensional time-domain wind speed sequence, where f i,t This represents the wind speed value at the i-th spatial point at the t-th time point; Angular frequency increment, Δω=ω up / N, where ω up The upper limit frequency; Re{·}: The operation of taking the real part of a complex number; The complex matrix after the inverse Fourier transform, and G in formula (10) (n×M) same; Time point index vector ; i: Imaginary unit, i 2=-1.

[0051] If the present invention adopts the JAX computing framework, it performs JIT just-in-time compilation optimization before reading parameters: sets the static_argnums static parameter for the single-frequency amplitude calculation function and the vectorized mapping operator, compiles it into machine code, and avoids the overhead of repeated compilation.

[0052] Step S2 represents the first level of parallel computation in this invention. Traditional methods, as shown in Algorithm 1, require nested loops to calculate the distance between all pairs of points in space. This invention eliminates nested loops through tensor broadcasting, reshaping the coordinate vectors into different dimensions (n×1 and 1×n). Utilizing the tensor broadcasting rules of modern numerical computing libraries (such as NumPy, JAX, and PyTorch), the vectors are automatically expanded into an n×n dimensional matrix during element-wise operations. A single vectorized operation can simultaneously compute all n... 2 A distance value that fully utilizes the CPU's SIMD instructions or the GPU's massively parallel thread capabilities.

[0053] Algorithm 1 Traditional Method: Calculating the distance matrix using nested loops 1: for i = 1 to n do 2: for j = 1 to n do 3:

[0054] 4: end for 5: end for Algorithm 2 Traditional Method: Sequential Loop Processing Frequency 1: for k = 1 to N do 2: Calculate the cross-spectral density matrix S k 3: Yes S k Perform Cholesky decomposition 4: Calculate the amplitude vector B k 5: end for

[0055] Steps S3 to S5 represent the second level of parallel computation in this invention. Traditional methods, as shown in Algorithm 2, require sequential loop processing of each frequency. This invention achieves the following key optimizations by defining a single-frequency amplitude calculation function, single_freq_amplitude: (1) The function uses tensor reshaping and tensor broadcasting mechanisms to achieve automatic dimension matching between spatial points, avoid loop traversal, and improve the spatial parallel efficiency of single frequency amplitude calculation. (2) The function only returns the final amplitude vector. Intermediate variables such as cross-spectral density matrix and Cholesky decomposition are released immediately after the calculation is completed inside the function, which effectively reduces memory usage. (3) Use the vectorized mapping operator vmap to perform frequency vectorized mapping on the single_freq_amplitude function. Each batch of frequency amplitude matrix calculations is completed, and the complete amplitude matrix is ​​formed by splicing them together after the batch calculations are completed. To adapt to the memory limitations of different computing devices, step S4 introduces an adaptive memory management mechanism. By estimating the memory requirements for frequency calculations in a single batch and comparing them with available memory, the optimal batch size is automatically determined. Each batch employs spatial-frequency dual-layer parallel computation, with no data dependency between batches, allowing for serial processing. After each batch is processed, the memory occupied by intermediate variables such as the cross-spectral density matrix is ​​released, ensuring the stable operation of the method under conditions of large spatial point counts and frequency bands.

[0056] This invention also provides a dual-layer parallel stochastic wind field rapid simulation system, comprising a set of computer program modules stored in the memory of an electronic device, executed by a processor calling hardware resources, supporting dynamic scheduling of CPU / GPU heterogeneous computing resources, supporting both GUI interaction and API calls as input methods, and supporting the export of results in multiple formats such as MAT, CSV, and CAD; used to implement the above-mentioned dual-layer parallel stochastic wind field rapid simulation method, including: 1. Parameter Input Module: Corresponding to step S1, it supports two input methods: GUI interactive interface and RESTful API. It can directly import CAD format spatial coordinate files, meteorological station measured wind parameters, and operating condition configuration files, automatically complete parameter format verification, and support batch import of multiple operating conditions.

[0057] 2. JIT compilation module: Enabled only under JAX backend. It automatically sets static parameters for single-frequency amplitude calculation function and vectorized mapping operator before the parameters are loaded, and compiles them into optimized machine code to avoid the overhead of repeated compilation at runtime.

[0058] 3. Spatial distance matrix calculation module: corresponding to step S2, calling the tensor operation interface of JAX / PyTorch to realize tensor reshaping and dimension broadcasting, and obtaining the coherent weighted distance matrix of all spatial point pairs in a single calculation.

[0059] 4. Single-frequency amplitude calculation module: corresponding to step S3 of the method, it encapsulates a pure function with no side effects, and all internal operations use tensor broadcasting mechanism to achieve parallelism in spatial dimensions.

[0060] 5. Parameter verification and matrix correction module: Embedded single-frequency amplitude calculation process, performs semi-positive definite correction on cross-spectral density matrix, and automatically triggers parameter verification if decomposition fails, sequentially checking the validity of average wind speed, spatial coordinates and coherence coefficient, and returning the index of abnormal parameters and modification prompts.

[0061] 6. Adaptive memory management module: corresponding to step S4, including memory reading unit, memory estimation unit and adaptive batching unit. The memory reading unit reads the available memory of the current device in real time through the operating system API. The memory estimation unit estimates the memory requirement of a single batch according to formula (5). The adaptive batching unit calculates the optimal batch size and completes the batching according to formula (6).

[0062] 7. Frequency Parallel Amplitude Matrix Calculation Module: Corresponding to step S5, the vectorized mapping operator is called to perform single-frequency amplitude calculations in batches, and batch splicing is automatically completed.

[0063] 8. Time-domain wind speed calculation module: corresponding to step S6, batch execution of inverse fast Fourier transform and phase modulation, supports export of results in multiple formats, built-in automatic verification function of spectral characteristics and correlation, and outputs standardized calculation report.

[0064] Example 1:

[0065] Example 1 provides a fast simulation method for stochastic wind fields based on a two-layer parallel architecture. It is implemented using the aforementioned two-layer parallel stochastic wind field fast simulation system and is used to verify the computational accuracy and statistical characteristics of the method. The wind speed parameters of a mountainous environment are used, specifically as follows: Test environment: Hardware: Intel i9-10900K CPU (supports AVX-512 SIMD instruction set), NVIDIA RTX 3090 24GB GPU (10496 CUDA cores), 32GB DDR4 memory; Software: JAX 0.4.20, CUDA 11.7; the traditional method benchmark uses the open-source wind field simulation library pycwind 1.2.0; The simulation parameters include: Space configuration: Number of spatial points: n=121; Spatial distribution: Uniformly distributed along the vertical direction, with a height range of z∈[80,320] meters; Average wind speed: calculated according to power law distribution. The wind speed at a reference height of 10m Roughness index ; Time and frequency configuration: Number of frequency bands: N=3000; Simulation duration: T = 600 seconds; Upper limit angular frequency: ω up =10πrad / s; Frequency increment: Δω = π / 300; Wind spectrum model: The wind spectrum adopts the Davenport downwind wind spectrum. The wind spectrum model is shown in formula (12). The autospectral density of each spatial point is calculated. Coherence function attenuation coefficient: C x =16,C y =6,C z =10.

[0066] Execution process: Perform step S1 of the method of this invention to obtain simulation parameters. Generate a spatial coordinate matrix P∈ 121×3 The x and y coordinates are both 0, and the z coordinate is linearly distributed within the range of 80 to 320 meters. Calculate the average wind speed vector U∈ at each point. 121 Based on the number of frequency bands N=3000 and the upper limit frequency ω up =10π, generating the angular frequency vector ω∈ 3000 angular frequency increment Generate a random phase matrix Φ∈ 3000×121 Each element is randomly sampled from a uniform distribution in [0, 2π]. After receiving the parameters, the processor automatically verifies their validity and confirms that there are no abnormalities. In step S2 of the method of this invention, the spatial distance matrix calculation module performs parallel calculation of the spatial coherent weighted distance matrix according to formula (1). Through tensor reshaping and tensor broadcasting mechanisms, the spatial coherent weighted distance matrix D∈ is calculated. 121×121 .

[0067] In step S3 of the method of this invention, a single-frequency amplitude calculation function, single_freq_amplitude, is defined. All calculations within the function are spatially parallelized using a tensor broadcasting mechanism, and the output is an amplitude vector B at a certain frequency. 121 .

[0068] Execute step S4 of the method of this invention to perform adaptive memory management. The floating-point data type adopts 32-bit single precision, d float =4 bytes. Let the available memory M be... avail =8GB, security level =1.5. Initial setting frequency and batch size N batch=3000, the adaptive memory management module calculates the single batch memory requirement according to formula (5): GB Due to M need <M avail No need to adjust batch size, keep N batch =3000.

[0069] In step S5 of the method of this invention, the frequency parallel amplitude matrix calculation module completes the amplitude calculation for all frequencies according to formula (7), performs frequency vectorization mapping vmap on the single_freq_amplitude function, performs parallel calculations on all frequency batches, and obtains the complete amplitude matrix B∈ 3000×121 .

[0070] In step S6 of the method of the present invention, the time-domain wind speed calculation module generates a time-domain wind speed sequence according to formulas (9) to (11) and outputs a MAT format file; including: zero-filling, inverse fast Fourier transform and phase modulation of the frequency domain amplitude matrix to generate a time-domain wind speed sequence f∈ 121×6000 .

[0071] Figure 2 The simulated wind speed time series are presented for different spatial locations. Each series reflects the fluctuation of wind speed over time at a specific spatial location. The overall temporal evolution and variability are consistent, reflecting the spatial correlation structure of the turbulent wind field. With increasing altitude, the amplitude of wind speed fluctuations increases, consistent with the expected pattern of power-law wind speed profiles. To assess the stationarity of the simulated time series, the Augmented Dickey-Fuller (ADF) test is applied to each series. Figure 3 The p-values ​​for the ADF stationarity test of wind speed series at all spatial locations were all well below the 0.05 significance level (the largest p-value was approximately 6 × 10⁻⁶). -4 This indicates that the simulated wind speed time series is stationary, consistent with the assumptions of the spectral model used.

[0072] To verify the accuracy of the wind spectrum simulation, the following wind field quality verification method can be selected: For the generated time-domain wind speed sequence f... (n×M) The automated accuracy verification process is as follows: For each spatial point, the wind speed time series is estimated using the Welch method to calculate the measured simulated power spectral density. The default Welch method parameters are: window length is 1 / 8 of the number of time points M, adjacent window overlap is 50%, and a Hamming window is added to reduce spectral leakage. The measured simulated power spectrum is compared with the target Davenport theoretical spectrum (Formula (12)) to calculate the average fitting error. Figure 4 As shown in af, Figure 4The af in the formula contains six independent subgraphs, and the horizontal axis of each independent subgraph is the dimensionless frequency f, which corresponds to the dimensionless frequency parameter defined in formula (12). The vertical axis of each independent subplot is the frequency-weighted normalized downwind power spectral density, which is obtained by multiplying the output of formula (12) by the physical frequency n. f get; Figure 4 The blue simulated spectrum curves of each independent subplot in f are obtained by power spectrum estimation using the Welch method from the wind speed time series generated by formula (11), and the theoretical spectrum curves are directly calculated by formula (12). The fitting error between the two in the full frequency range is less than 0.08%, which verifies the accuracy of wind spectrum simulation.

[0073] To further verify the spatial correlation structure, the cross-correlation function between different spatial points was calculated. For example... Figure 5 As shown in ad, Figure 5 The 'ad' plot contains four independent subplots. The horizontal axis of each independent subplot represents the time delay in seconds, and the vertical axis represents the spatial cross-correlation coefficient. The simulated cross-correlation curve of the measured blue line in the figure is calculated from the wind speed time series of two spatial points generated by formula (11) using the cross-correlation function. The theoretical cross-correlation curve of the orange line is calculated from the Davenport coherence function (formula (2)) using inverse Fourier transform. The average error between the two is less than 0.12%, which verifies the accuracy of the spatial correlation simulation. To facilitate direct comparison, both the simulated and theoretical correlation functions are normalized using their respective maximum absolute values. Figure 5 The `ad` parameter in the figure illustrates the comparison of cross-correlation functions between selected pairs of spatial points. In each subplot, the correlation function peaks at zero time delay and decays rapidly with increasing time delay, consistent with the expected behavior of turbulent wind fields. These results confirm that the method of this invention can accurately reproduce the spectral characteristics and spatial correlation features of the target wind field model.

[0074] Example 2:

[0075] Example 2 is used to evaluate the computational efficiency of the method of the present invention. Using the same spatial and temporal configuration parameters as Example 1, the computation time performance under different computing backends (NumPy, JAX, and PyTorch) was tested. The specific test scheme is as follows: Test group 1: Test the effect of different spatial point numbers n on computation time, with fixed frequency segment number N=3000, and the spatial point number n varies from 10 to 1000; Test Group 2: Test the effect of different frequency bands N on computation time, with the number of spatial points n=100 and the number of frequency bands N varying from 1000 to 10000; All other parameters and execution procedures are the same as in Example 1. Test results are as follows: Figures 6a-6b As shown.

[0076] Figure 6a The computation time performance under different numbers of space points is presented. The results show that computation time increases monotonically with increasing space points. The NumPy backend shows the most significant increase in computation time, mainly due to the limitations of its CPU-based computing architecture. In contrast, the JAX and PyTorch backends maintain a computation time of approximately 0.6 seconds when the number of space points is small (n≤200), demonstrating excellent scalability. As the number of space points further increases, the computation time gradually rises, mainly due to the activation of adaptive frequency batching mechanisms to cope with increased memory requirements. Even with the maximum space configuration n=1000, the computation time of both GPU backends is still less than 4 seconds, demonstrating a significant efficiency advantage compared to the hours or even days required by traditional methods.

[0077] Figure 6b The computation time performance at different frequency bands is presented. The results show that computation time increases monotonically with the number of frequency bands. The NumPy backend shows the most significant increase in computation time, while the JAX and PyTorch backends show relatively slower increases. This is because the JAX and PyTorch backends utilize GPU acceleration, fully leveraging the performance of modern computing hardware to improve the computational efficiency of wind field simulation.

[0078] The experimental results above demonstrate that the method of this invention can significantly improve the computational efficiency of stochastic wind field simulation while maintaining complete physical accuracy, making it suitable for simulation tasks with a large number of spatial points and frequency bands. Furthermore, the method of this invention has good scalability and can flexibly adapt to computational tasks of different scales.

[0079] The above description of the embodiments is provided to enable those skilled in the art to understand and use the present invention. It will be apparent to those skilled in the art that various modifications can be easily made to these embodiments, and the general principles described herein can be applied to other embodiments without inventive effort. Therefore, the present invention is not limited to the above embodiments, and any improvements and modifications made by those skilled in the art based on the disclosure of the present invention without departing from the scope of the present invention are within the protection scope of the present invention.

Claims

1. A fast simulation method for stochastic wind fields based on a two-layer parallel architecture, implemented using a harmonic superposition-type stochastic process simulation framework, characterized in that... The dual-layer parallelism includes two layers: spatial dimension parallelism and frequency dimension parallelism, and includes the following steps: S1: Obtain simulation parameters, including: number of spatial points n, number of frequency bands N, number of time points M, and spatial coordinates. Average wind speed vector Angular frequency vector Random phase matrix Angular frequency increment Spatial coherence attenuation coefficient Surface friction wind speed u * Where the angular frequency ω is related to the physical frequency n f Satisfying ω=2πn f ; S2: Parallel computation of the spatial coherence weighted distance matrix in spatial dimensions: The spatial coordinate vectors x, y, z are reshaped into n×1 and 1×n tensors, respectively, denoted as... and Using the tensor broadcasting mechanism, the matrix is ​​automatically expanded to an n×n dimension through element-level operations, and the spatial coherence weighted distance matrix is ​​obtained through parallel computation. The calculation formula is as follows: Official (1) In formula (1), : an n×n dimensional spatial coherent weighted distance matrix, where D i,j This represents the distance between the i-th spatial point and the j-th spatial point; Spatial coherence attenuation coefficients in the x, y, and z directions are used to describe the degree of attenuation of wind speed correlation in different directions; : The x, y, z coordinates of the i-th spatial point, with a tensor dimension of n×1; : The x, y, z coordinates of the j-th spatial point, with a tensor dimension of 1×n; S3: Defining a single frequency amplitude vector B (n) the spatial parallel computing function single_freq_amplitude comprising the following substeps: S31: Transform the average wind speed vector U z Remodeled into n×1 and 1×n dimensional tensors, denoted as U zi and U zj ; S32: Calculate the coherence function matrix: Official (2) In formula (2), : an n×n dimensional Davenport coherence function matrix, where Γ i,j This represents the coherence coefficient between the i-th spatial point and the j-th spatial point; ω: The currently calculated angular frequency value, which is an element in the angular frequency vector; D: Spatial coherence weighted distance matrix, which is the same as D in formula (1). (n×n) same; U zi : The average wind speed at the i-th spatial point, with a tensor dimension of n×1; U zj : The average wind speed at the j-th spatial point, with a tensor dimension of 1×n; S33: Calculate the autospectral density of each spatial point based on the Davenport downwind wind spectrum model, and reshape it into n×1 and 1×n dimensional tensors, denoted as S. i and S j ; S34: Calculate the cross-spectral density matrix: Official (3) In formula (3), : an n×n dimensional cross-spectral density matrix, where S i,j This represents the cross-spectral density between the i-th spatial point and the j-th spatial point; : Davenport coherence function matrix, and Γ in formula (2) (n×n) same; ⊙: Element-wise multiplication, multiplying corresponding elements; S i : The autospectral density of the i-th spatial point, with tensor dimension n×1; : The autospectral density of the j-th spatial point, where the superscript T is the matrix transpose operator and the tensor dimension is 1×n; S35: Perform Cholesky decomposition on the cross-spectral density matrix to obtain the lower triangular matrix H. (n×n) So that S=HH T ; S36: Calculate the amplitude vector at the current frequency: ； In formula (4), : n-dimensional amplitude vector, where B i This represents the amplitude of the i-th spatial point at the current frequency; An n×n lower triangular matrix, obtained by Cholesky decomposition of the cross-spectral density matrix, satisfying S=HH T ; i: Imaginary unit, i 2 =-1; : an n-dimensional random phase vector, where This represents the random phase of the i-th spatial point at the current frequency, with a value range of [0, 2π]. S4: Adaptive memory management: Dynamically determines the frequency and batch size, including: First, estimate the memory requirements for frequency calculation in a single batch, and then calculate the optimal batch size N based on the available device memory. batch ; Then divide the N frequencies into One batch, This is for rounding up; S5: Parallel computation of the complete amplitude matrix in the frequency dimension B (N×n) The single_freq_amplitude function defined in step S3 is vectorized using the vectorized mapping operator vmap, with a batch size of N for each calculation. batch The frequency is calculated, and after K batch calculations are performed, the complete amplitude matrix is ​​obtained by splicing them together; S6: Generate time-domain wind speed sequence: Perform transpose, zero-filling, inverse fast Fourier transform, and phase modulation on the frequency domain amplitude matrix to generate the time-domain wind speed sequence f. (n×M) The number of time points M satisfies M≥2N to avoid frequency domain aliasing.

2. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, In step S35, a positive semidefinite correction is performed on the cross-spectral density matrix before Cholesky decomposition to avoid task interruption due to decomposition failure. The positive semidefinite correction rule is as follows: When an accuracy better than 0.1% is required, the disturbance coefficient σ is taken as 10. -8 For conventional engineering calculations, σ = 10. -7 When the allowable accuracy error is ≥1%, take σ=10. -6 The correction method involves adding a small perturbation term σI to the cross-spectral density matrix, where I is the identity matrix; if σ=10 is used... -6 If the Cholesky decomposition still fails, the input parameter validation process is automatically triggered, returning the location of the abnormal parameters and modification prompts.

3. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, Step S4 includes the following sub-steps: S41: Single Batch Memory Requirements Calculation: Official (5) In formula (5), : Memory required for a single batch frequency calculation, in GB; n: number of space points, 3n 2 The memory usage corresponds to three types of n×n-dimensional tensors: spatial coherence weighted distance matrix, coherence function matrix, and cross-spectral density matrix. 4n corresponds to four types of n-dimensional tensors: average wind speed vector, self-spectral density vector, random phase vector, and amplitude vector. The number of bytes for floating-point data types: 4 bytes for 32-bit single-precision floating-point numbers and 8 bytes for 64-bit double-precision floating-point numbers; : Frequency of processing per batch; S42: Set safety factor Determine M need · Is it less than or equal to the available memory M? avail ,in For safety factors, the value ranges from 1.2 to 1.5; S43: If the limit is exceeded, calculate the optimal batch size using the following formula: Official (6) In formula (6), Optimal batch size, i.e., the frequency of processing per batch; : Round down; Available memory size, in GB; Safety factor, ranging from 1.2 to 1.5, is used to reserve a certain amount of memory space to avoid memory overflow; S44: Divide the N frequencies into One batch, for use in subsequent calculations.

4. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, The calculation method for the complete amplitude matrix in step S5 includes first calculating... The amplitude matrix is ​​transposed and then concatenated to form the complete amplitude matrix. : The The formula for calculating the amplitude matrix is: Official (7) In formula (7), : The amplitude matrix is ​​dimensional, where B k,i This represents the amplitude of the i-th spatial point at the k-th frequency; Vectorized mapping operator, used to apply a single-frequency amplitude calculation function in batches to multiple frequencies; : Single-frequency amplitude calculation function, the input is frequency value and random phase vector, and the output is amplitude vector; :N batch A 3D angular frequency vector containing all frequency values ​​processed in the current batch; : A dimensional random phase matrix, where Φ k,i This represents the random phase of the i-th spatial point at the k-th frequency; Will 3D amplitude matrix The transpose operation is ; The formula for splicing the complete amplitude matrix is: Official (8) In formula (8), : an n×N dimensional complete amplitude matrix, containing the amplitude of all spatial points at all frequencies; The amplitude matrix of the kth batch, with dimensions n×N batch ; K: Total number of batches, K= ,in This is for rounding up.

5. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, Step S6 includes sub-steps of zero-filling, inverse fast Fourier transform, and phase modulation: S61: Zero-filling: Concatenate an n×(MN) dimensional zero matrix to obtain the zero-filled amplitude matrix: Official (9) In formula (9), The amplitude matrix after zero-filling has a dimension of n×M; : The original amplitude matrix, with dimensions n×N; A zero matrix with n rows (MN) and columns, used to extend the number of columns of the amplitude matrix from N to M; M: Number of time points, i.e., the length of the generated time-domain wind speed sequence; S62: Inverse Fast Fourier Transform, the calculation formula is: Official (10) In formula (10), The complex matrix after the inverse Fourier transform has dimensions n×M; IFFT(·): Inverse Fast Fourier Transform operator, used to convert frequency domain signals into time domain signals; The amplitude matrix after zero-filling is the same as B in formula (9). padded same; S63: Phase modulation, calculated using the following formula: Official (11) In formula (11), : n×M dimensional time-domain wind speed sequence, where f i,t This represents the wind speed value at the i-th spatial point at the t-th time point; Angular frequency increment, Δω=ω up / N, where ω up The upper limit frequency; Re{·}: The operation of taking the real part of a complex number; The complex matrix after the inverse Fourier transform, and G in formula (10) (n×M) same; Time point index vector ; i: Imaginary unit, i 2 =-1.

6. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, In step S33, the formula for calculating the Davenport downwind wind spectrum is: Official (12) In formula (12), : Power spectral density of wind speed in the downwind direction; : Surface friction wind speed, in m / s; n f Physical frequency, measured in Hz; f: dimensionless frequency , where U is the average wind speed at height z in m / s, and z is the height above the ground in m.

7. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, The implementation of the vectorized mapping operator vmap in step S5 includes: (a) JAX backend: Calls the jax.vmap function and combines it with jax.jit for just-in-time compilation optimization; (b) PyTorch backend: Calls the torch.vmap function and uses CUDA for GPU acceleration.

8. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, In step S5, the amplitude matrix calculation for each batch of frequencies is performed in parallel. There is no data dependency between batches, and a serial processing method is used. After each batch is processed, the memory occupied by intermediate variables, including the cross-spectral density matrix and the lower triangular matrix of Cholesky decomposition, is released, effectively reducing the overall memory usage.

9. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, The inverse fast Fourier transform described in step S6 adopts a fully parallel computing method, performing inverse fast Fourier transform operations simultaneously on the frequency domain sequence of n spatial points.

10. The fast simulation method for stochastic wind fields based on a two-layer parallel architecture according to claim 1, characterized in that, If a JAX computing backend is used, step S0 is also included: before executing step S1, just-in-time (JIT) optimization is performed on the single-frequency amplitude calculation function and the vectorized mapping operator vmap to compile the function into efficient machine code in order to improve the overall computing performance.

11. A fast simulation system for stochastic wind fields based on a two-layer parallel architecture, characterized in that, The system is a collection of computer program modules stored in a memory and executed by a processor, used to implement the rapid simulation method for stochastic wind fields as described in any one of claims 1-10, including: (a) Parameter input module: used to acquire and store parameters required for simulation, including spatial coordinates, wind speed, frequency, phase, and coherence coefficient; (b) Spatial distance matrix calculation module: used to reshape spatial coordinates into tensors, achieve dimension matching through tensor broadcasting mechanism, calculate spatial coherent weighted distance matrix in parallel, and output the distance matrix of all spatial point pairs in one go without nested loops; (c) Single-frequency amplitude calculation module: Defines the spatial parallel calculation function for single-frequency amplitude vector, including steps for coherence function matrix calculation, cross-spectral density matrix calculation, Cholesky decomposition and amplitude vector calculation. The entire process adopts tensor broadcasting mechanism to achieve spatial parallelism. (d) Parameter verification and matrix correction module: used to perform positive semidefinite correction on the cross-spectral density matrix, and to trigger the input parameter verification process when the decomposition fails, returning the location of abnormal parameters and modification prompts; (e) Adaptive memory management module: including memory reading unit, memory estimation unit and adaptive batching unit. The memory reading unit is used to obtain the available memory of the current device through the operating system API. The memory estimation unit is used to add a frame correction factor to estimate the memory requirement of a single batch. The adaptive batching unit is used to calculate the optimal batch size and process the batches when the memory requirement exceeds the limit. (f) Frequency parallel amplitude matrix calculation module: Using the vectorized mapping operator vmap, the amplitude matrix of a frequency batch is calculated each time, and finally spliced ​​to obtain the complete amplitude matrix; (g) Time-domain wind speed calculation module: used to perform zero-filling, inverse fast Fourier transform and phase modulation on the frequency domain amplitude matrix to generate a time-domain wind field sequence.

12. The fast simulation system for stochastic wind fields based on a two-layer parallel architecture according to claim 11, characterized in that, It also includes a JIT compilation module: used in the JAX backend to perform just-in-time compilation optimization by setting static parameters for the single-frequency amplitude calculation function and the vectorized mapping operator vmap.

13. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the fast simulation method for stochastic wind fields based on a two-layer parallel architecture as described in any one of claims 1 to 10.

14. A computer-readable storage medium, characterized in that, It stores computer program instructions, which, when executed by a processor, implement the fast simulation method for stochastic wind fields based on a two-layer parallel architecture as described in any one of claims 1 to 10.

Images