A high-precision PINN multi-field coupling simulation method and system

Abstract

The application discloses a high-precision PINN multi-field coupling simulation method and system, relates to the technical field of numerical simulation, and comprises the following steps: obtaining space-time coordinates and initial values of stratum parameters and carrying out dimensionless scale reduction to eliminate the order-of-magnitude imbalance of multi-physical fields; building a full-connection deep neural network and enabling double-precision floating-point operation on the full link to avoid the operation truncation error of small parameters; with the aid of automatic differentiation, constructing a THM coupling residual according to a porous elastic control equation, and combining geometric feature perception to encrypt sampling of boundaries and corner points to suppress boundary numerical oscillation; setting asymmetric weights to construct a multi-objective loss function, and adopting Adam in combination with L-BFGS two-stage collaborative optimization to iteratively complete stratum permeability inversion. The application can improve simulation stability and parameter identification precision, and is suitable for multi-field simulation calculation of deep earth energy, carbon sequestration and other projects.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology, and more specifically to a high-precision PINN multi-field coupling simulation method and system. Background Technology

[0002] Thermo-hydro-mechanical (THM) coupling processes in porous media are widely present in critical engineering projects such as deep-earth energy extraction, geological disposal of nuclear waste, and carbon sequestration. Numerical simulation of these processes is of great significance for engineering safety assessment. Currently, although Physical Information Neural Networks (PINNs) have made some progress in fields such as fluid dynamics, they still face significant challenges in handling strongly coupled THM processes in porous media. (1) Multi-field gradient competition and magnitude imbalance: The fluid pressure field and the solid displacement field are deeply coupled through the effective stress principle. The physical units and magnitudes of the two are very different. Conventional PINN is prone to serious gradient competition during training, which leads to the loss function failing to converge.

[0003] (2) Numerical loss due to extremely small physical parameters: Core physical parameters such as formation permeability are usually as low as Order of magnitude. In a conventional 32-bit single-precision (float32) computing environment, such extremely small values ​​and their partial derivatives are highly susceptible to numerical truncation errors during backpropagation, causing the parameter inversion to deviate completely from the true value.

[0004] (3) Boundary layer numerical oscillation: Under complex 2D / 3D boundary conditions, the physical field changes drastically at the boundary. Usually, sampling strategies are difficult to capture local high gradient features, which often leads to non-physical oscillations in the simulation results near the boundary.

[0005] Therefore, how to solve the problems of inaccurate parameter inversion and non-convergence of calculation in multidimensional strongly coupled geological models is a technical problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0006] In view of this, the present invention provides a high-precision PINN multi-field coupling simulation method and system, solving the problems existing in the background technology. Specifically, the original spatial coordinates of the two-dimensional computational domain... Horizontal displacement of solids Vertical displacement of solid The physical magnitude of pore pressure differs greatly from that of pore pressure, easily leading to gradient competition and training non-convergence; abrupt changes in the physical field gradient at the upper and lower boundaries corresponding to model height and depth can easily cause boundary numerical oscillations; directly using the original coordinates for single-precision calculations can introduce truncation errors for minimal formation parameters, ultimately causing the inversion of key parameters such as permeability to deviate from the true value. To address these three technical pain points, this invention proposes a complete improvement scheme, including dimensionless scaling, end-to-end double-precision computation, geometrically aware non-uniform encrypted sampling, multi-objective asymmetric loss weighting, and two-stage collaborative optimization.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A high-precision PINN multi-field coupling simulation method includes the following steps: S1. Obtain the spatiotemporal coordinate data and the initial guess values ​​of the parameters to be inverted in the strata, and perform dimensionless scaling preprocessing on the spatiotemporal coordinate data and the physical output quantities to be solved. S2. Build a fully connected deep neural network and force the underlying operators and tensor data types to be declared as double-precision floating-point numbers. S3. Using automatic differentiation technology, the residual terms corresponding to the thermal-fluid-mechanical coupling law are constructed based on the porous elastic control equation; at the same time, non-uniform densification sampling is performed on the corner points and physical boundaries of the computational domain based on the geometric feature sensing method. S4. Using an asymmetric weighting configuration, construct a multi-objective integrated loss function that includes physical constraints, observation data, and boundary conditions; S5. Set the parameters of the formation to be identified as learnable parameters of the network, and use a two-stage optimization strategy to collaboratively iteratively update the network weights and the parameters to be identified until the preset convergence conditions are met. Then, output the multiphysics simulation results and the identified formation permeability parameters.

[0008] Optionally, in S1, the calculation formula for the dimensionless scaling preprocessing is:

[0009]

[0010]

[0011]

[0012] In the formula: For characteristic length, For characteristic time, For characteristic pressure, Characteristic displacement; For a two-dimensional computational domain Original spatial coordinates of direction for Dimensionless planar spatial coordinates of direction; , These represent the original time and the corresponding dimensionless quantity. , These are the formation fluid pressure and the corresponding dimensionless quantity, respectively. , These are the vertical displacement of the solid and the corresponding dimensionless quantities.

[0013] Optionally, in S2, the fully connected deep neural network adopts a decoupled multi-head network architecture, containing four independent sub-networks, corresponding to four types of physical field variables: pore pressure, horizontal displacement, vertical displacement, and volumetric strain. Each sub-network is coupled only at the loss function level through PDE residuals. The structure of a single sub-network includes an enhanced input layer, three hidden layers, and a single-neuron output layer. All network layers use the tanh activation function. In the original spatiotemporal coordinates... Building upon this foundation, the input layer is enhanced by injecting additional first-order frequency features, second-order frequency features, and decaying memory features, expanding the input dimension to 7-10 dimensions; among which, For a two-dimensional computational domain Original spatial coordinates of direction For a two-dimensional computational domain Original spatial coordinates of direction This is the original time.

[0014] Optionally, in S3, the residual terms corresponding to the thermo-fluid-mechanical coupling law include fluid diffusion residuals and solid equilibrium residuals, expressed as follows:

[0015]

[0016] In the formula: For fluid diffusion residuals, For solid equilibrium residuals; The pressure is a dimensionless formation fluid pressure. For dimensionless time, for Dimensionless planar spatial coordinates of direction for Dimensionless planar spatial coordinates of direction This represents the dimensionless vertical displacement of a solid. , where is the Biot coefficient, characterizing the contribution ratio of pore pressure to effective stress; This represents the ratio of material constants.

[0017] Optional, in S4, the comprehensive loss function The expression is:

[0018] In the formula: This is the physical residual loss term. The loss term is used to fit the observed data. The loss term is defined by boundary condition constraints. , , These are the weighting coefficients for the corresponding loss terms; The calculation methods for each loss item are as follows: The physical residual loss term uses the in-domain placement points to calculate the mean square error. The placement points adopt a hybrid sampling strategy, and the automatic differentiation engine is used to solve the partial derivatives of each order and calculate the average value of the residual sum of squares. The loss term for fitting the observation data is calculated based on the mean square error between the network prediction and the actual observation values ​​at the sensor observation points within the domain. The observation points are arranged along the cross profile at the center of the computational domain. The boundary condition constraint loss term consists of initial condition loss, Dirichlet boundary condition loss, and von Neumann and operator-type boundary condition loss. Meanwhile, a gradient pathology adaptive weighting mechanism is introduced to dynamically adjust the weight coefficients of each loss term to balance multi-objective optimization.

[0019] Optionally, in S5, the two-stage optimization strategy is as follows: First, a global search is performed using the Adam optimizer to initially determine the parameter range; then, the algorithm switches to the L-BFGS quasi-Newton algorithm for local high-precision fitting, and the penetration rate is monitored. Elastic modulus By observing the evolutionary trajectory of the iteration steps, the true values ​​of permeability and elastic modulus are finally obtained.

[0020] A high-precision PINN multi-field coupling simulation system, executing a high-precision PINN multi-field coupling simulation method as described in any one of the above, including: The input module is used to acquire spatiotemporal coordinate data and initial guesses of the parameters to be inverted in the formation; The preprocessing module is used to perform dimensionless scaling preprocessing on the spatiotemporal coordinate data and the physical output quantities to be solved. The network configuration module is used to build a fully connected deep neural network and force the underlying operators and tensor data types to be declared as double-precision floating-point numbers. The residual construction module is used to construct the residual terms corresponding to the thermo-fluid-mechanical coupling law based on the porous elastic control equation using automatic differentiation technology. The sampling module is used to perform non-uniform encrypted sampling of the corner points and physical boundaries of the computational domain based on geometric feature perception. The loss function construction module is used to construct a multi-objective integrated loss function that includes physical constraints, observation data, and boundary conditions using an asymmetric weight configuration method. The parameter optimization module is used to set the parameters of the formation to be identified as learnable parameters of the network, and adopts a two-stage optimization strategy to collaboratively and iteratively update the network weights and the parameters to be identified. The output module is used to output the multiphysics simulation results and the identified formation permeability parameters after the iteration meets the preset convergence conditions.

[0021] As can be seen from the above technical solution, compared with the prior art, the present invention discloses a high-precision PINN multi-field coupling simulation method and system, which has the following beneficial effects: (1) The present invention adopts a dimensionless scaling strategy in the multidimensional spatiotemporal domain, and reconstructs 2D / 3D spatiotemporal coordinates and physical outputs through feature factors, mapping physical quantities of different magnitudes to a unified numerical range, thereby eliminating gradient imbalance caused by differences in magnitude.

[0022] (2) The present invention adopts a full-path double-precision (float64) tensor operation control mechanism. In the entire process of neural network initialization, forward inference and automatic differentiation, hardware-level double-precision floating-point operation is forced to accurately capture the weak gradient changes caused by extremely small parameters and eliminate truncation error.

[0023] (3) In response to the inverse problem, this invention introduces a weight offset mechanism (such as increasing the weight of the data item to more than 100) to resolve the gradient conflict between the physical residual item and the observed data item through asymmetric weighting, thereby guiding the model to converge quickly under sparse observation.

[0024] (4) This invention performs non-uniform encrypted sampling on the corner points and physical boundaries of the computational domain, introduces boundary reinforcement constraints, effectively suppresses the numerical instability of the boundary layer, and ensures the physical consistency of the inversion parameters at geometric singularities. Attached Figure Description

[0025] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0026] Figure 1 The overall flowchart of the high-precision PINN multi-field coupling simulation method for identifying extremely small permeability parameters provided by the present invention is shown below. Figure 2 The structural diagram of the high-precision PINN multi-field coupling simulation system for identifying extremely low permeability parameters provided by the present invention is shown. Figure 3Comparison curves of calculation results of multiphysics field distribution along depth at different times under two floating-point precisions, single-precision float32 and double-precision float64, provided by the present invention. Figure 4 The present invention provides a comparison curve of the convergence of the total training loss function under different floating-point precisions. Figure 5 The curve showing the convergence change of penetration rate with training iterations under the Barry-Mercer example provided by this invention; Figure 6 The spatiotemporal distribution cloud map of the physical field and absolute error of the 2D Barry-Mercer example provided by this invention. Detailed Implementation

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

[0028] This invention discloses a high-precision PINN multi-field coupling simulation method, such as... Figure 1 As shown, it includes the following steps: S1. Obtain the spatiotemporal coordinate data and the initial guess values ​​of the parameters to be inverted in the strata, and perform dimensionless scaling preprocessing on the spatiotemporal coordinate data and the physical output quantities to be solved. S2. Build a fully connected deep neural network (MLP) and force the underlying operators and tensor data types to be declared as double-precision floating-point numbers (float64). S3. Using automatic differentiation (AD) technology, the residual terms corresponding to the thermal-fluid-mechanical coupling law are constructed based on the porous elastic control equation. At the same time, based on the geometric feature sensing method, non-uniform and dense sampling is performed on the corner points and physical boundaries of the computational domain to correct the numerical oscillations at the boundaries of the high-dimensional space. S4. Using an asymmetric weighting configuration, construct a multi-objective integrated loss function that includes physical constraints, observation data, and boundary conditions; S5. Set the parameters of the formation to be identified as learnable parameters of the network, and use a two-stage optimization strategy to collaboratively iteratively update the network weights and the parameters to be identified until the preset convergence conditions are met. Then, output the multiphysics simulation results and the identified formation permeability parameters.

[0029] Furthermore, addressing the issue of gradient imbalance caused by the significant difference in magnitude (typically 3-6 orders of magnitude) between the pressure and displacement fields in porous media, S1 of this embodiment proposes dimensionless processing of the 2D / 3D spatiotemporal coordinates and physical output. The calculation formula for the dimensionless scaling preprocessing is as follows:

[0030]

[0031]

[0032]

[0033] In the formula: The feature length is typically taken as the maximum geometric dimension of the computational domain; The characteristic time is usually taken as the characteristic response time of the physical process; The characteristic pressure is usually taken as the maximum value of the external load, the initial pore pressure difference, or the typical pressure level of the project. The characteristic displacement is usually taken as the maximum allowable deformation in the computational domain or the theoretical maximum displacement; For a two-dimensional computational domain Original spatial coordinates of direction for Dimensionless planar spatial coordinates of direction; , These represent the original time and the corresponding dimensionless quantity. , These are the formation fluid pressure and the corresponding dimensionless quantity, respectively. , These are the vertical displacement of the solid and the corresponding dimensionless quantities. This step transforms the original dimensional parameters into dimensionless parameters through scaling, unifying the numerical ranges of various physical fields, and fundamentally solving the gradient competition and model training non-convergence problems caused by the imbalance between the magnitudes of the pressure field and the displacement field.

[0034] Through the above transformation, a feature length is introduced. Characteristic Time Characteristic pressure and characteristic displacement Mapping the original physical field to a dimensionless space The aim is to eliminate the gradient vanishing problem in multi-field coupling by reconstructing the input space, thereby eliminating the difference of up to 3-6 orders of magnitude between the pressure field and the displacement field.

[0035] For penetration rates as low as 10 -14The extremely small numerical values ​​mean that conventional 32-bit single-precision numbers cannot provide enough significant bits, leading to severe truncation errors during backpropagation. Therefore, S2 in this embodiment proposes to force the underlying operators and tensor data types to be declared as double-precision floating-point numbers, ensuring accurate capture of weak physical gradients through hardware-level precision control.

[0036] In S2, the fully connected deep neural network adopts a decoupled multi-head network architecture, containing four independent subnetworks, corresponding to four types of physical field variables: pore pressure, horizontal displacement, vertical displacement, and volumetric strain. Each subnetwork is coupled only at the loss function level via PDE residuals. The structure of a single subnetwork includes an enhanced input layer, three hidden layers, and a single-neuron output layer. All network layers use the tanh activation function. (In the original spatiotemporal coordinates...) Building upon this foundation, the input layer is enhanced by injecting additional first-order frequency features, second-order frequency features, and decaying memory features, expanding the input dimension to 7-10 dimensions; among which, For a two-dimensional computational domain Original spatial coordinates of direction For a two-dimensional computational domain Original spatial coordinates of direction This is the original time.

[0037] Specifically, the fully connected deep neural network used in this embodiment has made clear improvements and optimizations compared to existing technologies: (1) Infrastructure (1D Tersaqi problem) The design of a fully connected feedforward neural network (MLP) is as follows: ① Input layer: 2 neurons, which receive dimensionless spatial coordinates and dimensionless time coordinates respectively; ② Hidden layers: 4 layers, 32 neurons per layer; Experimental comparison shows that compared with 20 neurons, the use of 32 neurons in this embodiment significantly improves both convergence speed and final accuracy; ③ Output layer: 2 neurons, which output dimensionless pore pressure and dimensionless displacement respectively; ④ Activation function: tanh (hyperbolic tangent function), used uniformly in all hidden layers; ⑤ Weight initialization: Glorot normal distribution initialization.

[0038] Overall structure: [2]→

[32] →

[32] →

[32] →

[32] →[2].

[0039] (2) Innovative architecture (2D Barry-Mercer problem) -- Decoupling multi-head network

[0040] To address the gradient competition problem in 2D multiphysics coupling problems, this embodiment proposes a decoupled multi-head network architecture. A significant improvement over existing technologies (such as standard single fully connected networks) is that it integrates different physical field variables (formation fluid pressure...) into a single network. Horizontal displacement of solids Vertical displacement of solid Volumetric strain ) allocated to 4 independent subnets (net_p, net_u, net_v, net_ Each network has its own independent network parameters, and they are coupled only at the loss function level through PDE residuals.

[0041] The structure of each sub-network is: input layer → hidden layer (128 neurons) → hidden layer (128 neurons) → hidden layer (128 neurons) → output layer (1 neuron), with the activation function being tanh.

[0042] Input feature enhancement specifically refers to: in standard spatiotemporal coordinates Based on this, time-phase features are injected, including: first-order frequency features. , Second-order frequency characteristics , Decaying memory characteristics , , This expands the 3D input to a 7-10 dimensional augmented input space. Among these, The angular frequency is a characteristic of the time phase, corresponding to the time periodicity of the physical process.

[0043] The advantages of the fully connected deep neural network in this embodiment compared to existing technologies include the following aspects: ① The decoupled architecture effectively alleviates the gradient competition problem between multiple physics fields—different physical quantities no longer interfere with each other by sharing network parameters; ② The injection of time phase features makes it easier for the network to learn periodic / decaying time responses, solving the problem that standard MLPs have difficulty learning high-frequency / periodic time components; ③ Compared to a single network architecture, decoupled networks compute each physical quantity independently during forward propagation, and only couple through PDE during backward propagation differentiation. This is equivalent to a soft operator splitting strategy, which effectively improves training stability.

[0044] Furthermore, in S3, the residual terms corresponding to the thermo-fluid-mechanical coupling law include fluid diffusion residuals and solid equilibrium residuals, expressed as:

[0045]

[0046] In the formula: For fluid diffusion residuals, For solid equilibrium residuals; The pressure is a dimensionless formation fluid pressure. For dimensionless time, for Dimensionless planar spatial coordinates of direction for Dimensionless planar spatial coordinates of direction This represents the dimensionless vertical displacement of a solid. , where is the Biot coefficient, characterizing the contribution ratio of pore pressure to effective stress; This represents the material constant ratio, specifically a combination of Poisson's ratio parameters, used in the dimensionless solid equilibrium equations. This step transforms physical laws into loss terms for the neural network, enabling the model to strictly adhere to the laws of conservation of mass and momentum even without massive amounts of labeled data.

[0047] In S3, "geometric feature perception" is one of the key technical features of this invention, and its determination principle and specific implementation method are as follows: (I) Geometric Singularity Identification and Avoidance Strategies Geometric singularities (such as corners and intersections of boundary lines) in the computational domain are subject to multiple boundary conditions simultaneously, leading to non-uniqueness or strong discontinuity in the physical field at these locations. This embodiment employs the following strategy to identify and avoid singularities: 1. Space-Time Buffer Zone Identification principle: In the spatiotemporal domain, when a point simultaneously satisfies two or more boundary conditions, that point is identified as a corner singularity (e.g., ...). (It is simultaneously constrained by the top surface boundary conditions and initial conditions).

[0048] Solution: Introduce a safe distance during sampling. ( (for time offset) and ( (For spatial offset), so that the boundary sampling points and the initial sampling points maintain a safe distance in the spatiotemporal domain, avoiding condition conflicts at singularities.

[0049] 2. Boundary Layered Sampling Strategy

[0050] Each boundary is sampled independently, with each having its own independent set of sampling points. Time sampling of the top boundary is avoided. (Through epsilon offset), initial conditional space sampling avoids (via epsilon offset), bottom boundary and There are no singularity collisions, so complete sampling is maintained. For the 1D Terzaghi problem, (bottom), (Top face); For the 2D Barry-Mercer problem, , ;in, , The computational domain is respectively in Dimensionless minimum and maximum coordinate values ​​of the direction (depth direction). , The computational domain is respectively in The dimensionless minimum and maximum coordinate values ​​of the direction.

[0051] (ii) Non-uniform dense sampling of the boundary layer

[0052] The physical field changes drastically near the boundary, making it difficult for uniform sampling to capture local high-gradient features, which can easily lead to numerical oscillations. The encryption strategy in this embodiment includes: 1. Source-Concentrated Gaussian Sampling Identification principle: Locate the planar coordinates of the injection source / load application point within the dimensionless computational domain. For example, in the classic Barry-Mercer example, the injection source is located in the lower left corner of the computation domain. .

[0053] Processing method: In the sampling within the domain, 50% of the configuration points are sampled using a Gaussian distribution with the source point as the center and sigma=0.04 as the variance, so that the high gradient region (near the source) can obtain a higher sampling density.

[0054] 2. Time-Phase Balanced Sampling

[0055] Identification principle: Divide the time domain evenly into 8 phase intervals (bins) to ensure that each phase interval has representative sampling points.

[0056] Processing method: Bucket equalization sampling is used for both boundary conditions and points within the domain to avoid training focusing on a certain time phase while ignoring other phases (which is especially important in periodic load problems).

[0057] 3. Residual-based Adaptive Refinement (RAR)

[0058] Identification principle: During the training process, a large number of candidate points (e.g., 2000) are randomly sampled periodically (e.g., every 4000 iterations) across the entire domain, and the absolute value of the PDE residual is calculated.

[0059] Processing method: Select the top one with the largest residual. Adding 200 points to the training set enables adaptive encryption in high residual regions. For the inverse problem, phase-bucketed RAR is further employed to ensure that the quota of newly added points is equally allocated to each time phase interval.

[0060] This step involves non-uniformly refined sampling of the geometric boundaries and corners corresponding to the model's height and depth. This accurately captures the high gradient features of the boundaries, eliminates non-physical oscillations in the boundary layer caused by uniform sampling of the original coordinates, and avoids distortion in the inversion of formation parameters.

[0061] Furthermore, in S4, the comprehensive loss function The expression is:

[0062] In the formula: This is the physical residual loss term. The loss term is used to fit the observed data. The loss term is defined by boundary condition constraints. , , These are the weighting coefficients for the corresponding loss terms.

[0063] The calculation methods for each loss item are as follows: ① The physical residual loss term uses the mean square error calculated by the in-domain placement points. The placement points adopt a hybrid sampling strategy, and the partial derivatives of each order are solved by the automatic differentiation engine to calculate the average value of the residual sum of squares.

[0064] Specifically, random sampling within the computational domain For each collocation point, the residuals of each PDE are calculated using Automatic Differentiation (AD) technology, and then the mean squared error (MSE) is calculated. The expression is:

[0065] In the formula: To calculate the total number of configuration points within the domain, the value is between 10,000 and 20,000. For the first The spatiotemporal coordinates of each sampling point For the fluid diffusion equation residuals, for Directional solid mechanics equilibrium equation residuals for Directional solid mechanics equilibrium equation residuals The residuals of the volumetric strain correlation equations, For volumetric strain.

[0066] The sampling of configuration points adopts a hybrid strategy: 40% uniform sampling + 20% Gaussian concentrated sampling of source points + 40% time phase equalization sampling.

[0067] Automatic differentiation is implemented using PyTorch / TensorFlow's autograd engine, which can accurately calculate partial derivatives of any order.

[0068] For the inverse problem, Residual Adaptive Encryption (RAR) is enabled during training to periodically supplement high residual points.

[0069] ② The loss term for fitting the observation data is calculated based on the mean square error between the network predicted value and the actual observed value, with the observation points arranged along the cross profile at the center of the computational domain.

[0070] Specifically, at the sensor / observation point location, the mean square error between the network's predicted value and the known observed value is calculated, expressed as:

[0071] In the formula: The total number of observation points. For the first Pore ​​pressure predicted by the network at each observation point. For the first The actual observed / measured pore pressure at each observation point. , For the horizontal and vertical displacements predicted by the network, , The measured horizontal and vertical displacements at the corresponding points.

[0072] Observation point arrangement: within the computational domain, along (Centerline) and Discrete sensors are arranged along the cross-section of the (centerline), with typical locations as follows: .in, To standardize the dimensionless reference length, used to scale spatial coordinates to... Dimensionless interval. For the 1D Terzaghi consolidation problem. Take the depth of the computational domain For the 2D Barry-Mercer problem, Take the side length of the computational domain (at this time (The computational domain is a square).

[0073] The observation data can be actual engineering monitoring data, calculation results of high-fidelity numerical simulation (such as the finite element method), or analytical solutions, such as Terzaghi one-dimensional consolidation theory solutions and Barry-Mercer series solutions.

[0074] For the inverse problem, the weight of the observed data is increased to 100 or more (asymmetric weighting strategy) to guide the model to quickly approximate the true value under sparse observations.

[0075] ③ The boundary condition constraint loss term consists of initial condition loss, Dirichlet boundary condition loss, and Neumann and operator-type boundary condition loss.

[0076] Specifically, sampling is performed at the spatial boundary and initial time surface of the computational domain, and the mean square error of the boundary condition residuals is calculated, expressed as:

[0077] In the formula: Loss due to initial conditions Loss at the Dirichlet border, This is a Neumann / operator-type boundary loss.

[0078] A. Initial conditions ( ): (Initial pore pressure is zero); (Initial displacement is zero); Number of sampling points: .

[0079] B. Dirichlet boundary conditions (specified value conditions): (All outer boundaries, drainage conditions); (Bottom (constraining vertical displacement). ( and (place); Number of sampling points: per edge One point.

[0080] C. Neumann / Operator boundary conditions (specifying gradient conditions): (Bottom (No drainage conditions) ( and (no horizontal strain) ( and (No vertical strain) Top surface stress boundary conditions: External load constraints are applied through operator form BC; Implementation: Use dde.icbc.OperatorBC to calculate gradients and apply constraints through automatic differentiation.

[0081] ④ Gradient Pathology (GP)

[0082] In addition to the three basic loss terms mentioned above, this embodiment introduces a gradient pathology adaptive weighting mechanism (GradientBalance / GradNorm) to dynamically adjust the weight coefficients of each loss term in order to balance multi-objective optimization.

[0083] The calculation steps are as follows: Every 100 iterations, 2048 points are randomly sampled from the training set; Calculate the gradient norm of each PDE residual with respect to the weights of each sub-network output layer; Using the median of the gradient norm as a reference: ;in, To smooth out the exponent and mitigate drastic changes in weights; Indicates the first The gradient norm corresponding to the residual loss term of the PDE class. The median of the set of all PDE gradient norms. For the first The target corrected weights for each PDE loss term; The target weights are pruned to prevent weight divergence; the inverse problem is limited to between 0.5 and 20, while the forward problem can be relaxed to between 1 and 1000. Update weights using an exponential moving average: ;in, This represents the original weights of the loss term used in the previous iteration. This indicates the ideal weight that is expected to be adjusted in this round. This represents the final weight after the exponential moving average update; Only the weights of the PDE residuals are adjusted, while the weights of the boundary conditions remain fixed.

[0084] This mechanism ensures that the gradient contributions of each PDE equation remain on the same order of magnitude, avoiding the problem of one equation (such as the fluid diffusion equation) dominating training while another equation (such as the momentum conservation equation) is ignored (i.e., the "silence" problem).

[0085] The core innovation of this invention in solving the problem of "inaccurate parameter inversion" is the construction of a hybrid loss function with asymmetric weights in S4. In particular, in response to the convergence stagnation caused by the dominance of the physics term in the inversion problem, this embodiment increases the weight of the observation data term to 100 or more.

[0086] Furthermore, in S5, the two-stage optimization strategy is as follows: First, a global search is performed using the Adam optimizer to initially determine the parameter range; then, the algorithm switches to the L-BFGS quasi-Newton algorithm for local high-precision fitting, and the penetration rate is monitored. Elastic modulus By observing the evolutionary trajectory of the iteration steps, the true values ​​of permeability and elastic modulus are finally obtained.

[0087] In practice, this embodiment can also employ a more refined "staged sequential training strategy" than the "two-stage" approach, with the following stage division and switching conditions: (a) Three-stage sequential training strategy (for the inverse problem) Stage 0 – Field Variable Warm-up: Number of iterations: 2000; Optimizer: Adam, learning rate 1e-3; Parameter status to be identified: Frozen ( (Not participating in training) Objective: To train only the weights of a neural network so that the network can quickly establish a reasonable physical field distribution; Loss weight configuration: PDE item [1,1,1,1], BC item

[50] ×12, data item

[50] ×6; Switching condition: Automatically enter the next stage after completing 2000 iterations.

[0088] Phase 1 - Strong Stage 1: Strong Identification Number of iterations: 6000; Optimizer: Adam, learning rate 5e-4; Status of parameter to be identified: thawed ( (Participate in training) Objective: To initiate the identification of permeability parameters based on the initial convergence of the physical field; Loss weighting configuration: PDE term [20,1,1,1] (The weight of the fluid mass conservation term is increased by 20 times to highlight permeability) The role of), BC item

[50] ×12, data item [150,150,20,20,20,20] (pressure observation data weight increased to 150); Switching condition: Automatically enter the next stage after completing 6000 iterations.

[0089] Stage 2 – Fully Coupled Joint Training: Number of iterations: 4000; Optimizer: Adam, learning rate 2e-4; Parameter status to be identified: Continue training; Objective: To fine-tune network weights and parameters to be identified in a coordinated manner to approximate the global optimum; Loss weight configuration: PDE item [5,5,5,2], BC item

[50] ×12, data item

[80] ×6; GP gradient pathology adaptive weights: Enabled (only in this phase), dynamically balancing each PDE term; Switching condition: Enter the L-BFGS fine-tuning stage after completing 4000 iterations.

[0090] (ii) Adam→L-BFGS switching

[0091] After Adam completes all training phases (approximately 12,000 iterations), it automatically switches to the L-BFGS quasi-Newton algorithm for local high-precision fitting. The switching criterion is: Rigidity determination: Complete the preset total number of Adam iterations (Stage0 + Stage1 + Stage2); Auxiliary monitoring indicators (used to determine whether to switch to Adam early or increase Adam iterations): (1) The loss function decreases gradually (relative decrease of <1% after 1000 consecutive iterations); (2) Parameters to be identified The evolutionary trajectory tends to stabilize (change <0.1% in 500 consecutive iterations); (3) The relative magnitudes of the PDE residuals no longer change significantly.

[0092] Reason for switching to L-BFGS: (1) L-BFGS does not require manual setting of the learning rate; the step size is automatically determined by the Hessian matrix. (2) L-BFGS has a superlinear convergence rate when approaching a local optimum; (3) Adam is responsible for global search (escaping saddle points and local extrema), and L-BFGS is responsible for local high-precision fitting.

[0093] Note: For positive problems (with known parameters), it is simplified to a direct two-stage switching from Adam to L-BFGS, with Adam switching after approximately 15,000 iterations.

[0094] and Figure 1 Corresponding to the method described above, this embodiment of the invention also provides a high-precision PINN multi-field coupling simulation system for... Figure 1 The specific implementation of the method, the high-precision PINN multi-field coupling simulation system provided in this embodiment of the invention, can be applied to computer terminals or various mobile devices, such as... Figure 2 As shown, it specifically includes: The input module is used to acquire spatiotemporal coordinate data and initial guesses of the parameters to be inverted in the formation; The preprocessing module is used to perform dimensionless scaling preprocessing on the spatiotemporal coordinate data and the physical output quantities to be solved. The network configuration module is used to build a fully connected deep neural network and force the underlying operators and tensor data types to be declared as double-precision floating-point numbers. The residual construction module is used to construct the residual terms corresponding to the thermo-fluid-mechanical coupling law based on the porous elastic control equation using automatic differentiation technology. The sampling module is used to perform non-uniform encrypted sampling of the corner points and physical boundaries of the computational domain based on geometric feature perception. The loss function construction module is used to construct a multi-objective integrated loss function that includes physical constraints, observation data, and boundary conditions using an asymmetric weight configuration method. The parameter optimization module is used to set the parameters of the formation to be identified as learnable parameters of the network, and adopts a two-stage optimization strategy to collaboratively and iteratively update the network weights and the parameters to be identified. The output module is used to output the multiphysics simulation results and the identified formation permeability parameters after the iteration meets the preset convergence conditions.

[0095] Figure 3 This demonstrates the formation fluid pressure in a 1D Terzaghi consolidation problem after training with float32 and float64 accuracies, respectively. (pore pressure) and solid vertical displacement The curve comparing the predicted values ​​with the analytical solutions. Specific content includes: Figure 3 First row: Pore pressure (MPa) along depth The distribution of (m) at 3 typical times ( , , The pore pressure distribution curves along the depth are plotted respectively, with the analytical solution plotted as a solid line, the float32 single-precision prediction result plotted as a dashed line, and the float64 double-precision prediction result of this invention plotted as a dotted-dash line. The absolute error distribution between the predicted value and the analytical solution is shown for each curve. Figure 3 Second bottom row: Vertical displacement of solid (m) along depth Corresponding comparison of (m); Quantitative conclusions: Regarding Figure 3 As shown in the 1D Terzaghi solid-state calculation example, the full-link float64 double-precision scheme of this invention reduces the full-field relative L2 prediction error of pore pressure and displacement field by approximately 13% and 38% respectively compared to the conventional float32 single-precision scheme; at the same time, the final total training loss is reduced by approximately 40 times (from...). Down to This indicates that the float64 scheme has significant advantages in convergence depth and optimization quality.

[0096] Figure 4 The graph separately displays the evolution curves of the total training loss with the number of iterations for two calculation methods: float32 single precision and float64 double precision. Using a semi-logarithmic coordinate system (logarithmic on the vertical axis and linear on the horizontal axis), the graph allows for a direct comparison of the convergence speed and the difference in the steady-state minimum loss value between the two floating-point precision models. The curve comparison shows that, at the end of the training iteration, the steady-state total loss of the float64 scheme is only 1 / 10 to 1 / 100 of that of the float32 scheme.

[0097] Figure 5 This demonstrates the parameters to be identified in the 2D Barry-Mercer inverse problem. The evolution trajectory of (dimensionless penetration rate) with the number of training iterations. Specific details include: Solid line: PINN's prediction The evolution of the value with the number of iterations; dotted line: The true value (baseline); The vertical axis represents permeability. The physical value (m²) is expressed in scientific notation; The horizontal axis represents the number of training iterations; based on Figure 5 It can be observed that: Initial stage The guessed value (deviating from the true value) gradually converges to near the true value, and the relative error is eventually within an acceptable range (e.g., <5%).

[0098] Figure 6This paper presents a comparison of the physical field distribution near the boundary under different sampling strategies in the 2D Barry-Mercer problem. The figure includes two sets of two-dimensional contour maps for pore pressure and displacement. Each set has three subplots: the analytical solution, the PINN predicted solution, and the absolute error contour map between the predicted and analytical solutions. A contour map with color filling is used to display the two-dimensional full-field distribution, allowing for intuitive observation of the error distribution characteristics in the boundary region. After optimization using geometric feature-aware encrypted sampling, the prediction results in the boundary region no longer exhibit significant non-physical numerical oscillations.

[0099] In addition, specific observation test points within the computational domain can be selected to plot time series comparison curves of predicted values ​​and analytical solutions, further verifying the effect of the boundary sampling optimization strategy on improving the overall prediction accuracy.

[0100] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the systems disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the descriptions are relatively simple; relevant parts can be referred to the method section.

[0101] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A high-precision PINN multi-field coupling simulation method, characterized in that, Includes the following steps: S1. Obtain the spatiotemporal coordinate data and the initial guess values ​​of the parameters to be inverted in the strata, and perform dimensionless scaling preprocessing on the spatiotemporal coordinate data and the physical output quantities to be solved. S2. Build a fully connected deep neural network and force the underlying operators and tensor data types to be declared as double-precision floating-point numbers. S3. Using automatic differentiation technology, the residual terms corresponding to the thermal-fluid-mechanical coupling law are constructed based on the porous elastic control equation; at the same time, non-uniform densification sampling is performed on the corner points and physical boundaries of the computational domain based on the geometric feature sensing method. S4. Employ an asymmetric weighting approach to construct a multi-objective integrated loss function that incorporates physical constraints, observation data, and boundary conditions; integrated loss function The expression is: In the formula: This is the physical residual loss term. The loss term is used to fit the observed data. The loss term is defined by boundary condition constraints. , , These are the weighting coefficients for the corresponding loss terms; The calculation methods for each loss item are as follows: The physical residual loss term uses the in-domain placement points to calculate the mean square error. The placement points adopt a hybrid sampling strategy, and the automatic differentiation engine is used to solve the partial derivatives of each order and calculate the average value of the residual sum of squares. The loss term for fitting the observation data is calculated based on the mean square error between the network prediction and the actual observation values ​​at the sensor observation points within the domain. The observation points are arranged along the cross profile at the center of the computational domain. The boundary condition constraint loss term consists of initial condition loss, Dirichlet boundary condition loss, and von Neumann and operator-type boundary condition loss. Meanwhile, a gradient pathology adaptive weighting mechanism is introduced to dynamically adjust the weight coefficients of each loss term in order to balance multi-objective optimization. S5. Set the parameters of the formation to be identified as learnable parameters of the network, and use a two-stage optimization strategy to collaboratively iteratively update the network weights and the parameters to be identified until the preset convergence conditions are met. Then, output the multiphysics simulation results and the identified formation permeability parameters.

2. The high-precision PINN multi-field coupling simulation method according to claim 1, characterized in that, In S1, the calculation formula for the dimensionless scaling preprocessing is: In the formula: For characteristic length, For characteristic time, For characteristic pressure, Characteristic displacement; For a two-dimensional computational domain Original spatial coordinates of direction for Dimensionless planar spatial coordinates of direction; , These represent the original time and the corresponding dimensionless quantity. , These are the formation fluid pressure and the corresponding dimensionless quantity, respectively. , These are the vertical displacement of the solid and the corresponding dimensionless quantities.

3. The high-precision PINN multi-field coupling simulation method according to claim 1, characterized in that, In S2, the fully connected deep neural network adopts a decoupled multi-head network architecture, containing four independent subnetworks, corresponding to four types of physical field variables: pore pressure, horizontal displacement, vertical displacement, and volumetric strain. Each subnetwork is coupled only at the loss function level via PDE residuals. The structure of a single subnetwork includes an enhanced input layer, three hidden layers, and a single-neuron output layer. All network layers use the tanh activation function. (In the original spatiotemporal coordinates...) Building upon this foundation, the input layer is enhanced by injecting additional first-order frequency features, second-order frequency features, and decaying memory features, expanding the input dimension to 7-10 dimensions; among which, For a two-dimensional computational domain Original spatial coordinates of direction For a two-dimensional computational domain Original spatial coordinates of direction This is the original time.

4. The high-precision PINN multi-field coupling simulation method according to claim 1, characterized in that, In S3, the residual terms corresponding to the thermo-fluid-mechanical coupling law include fluid diffusion residuals and solid equilibrium residuals, expressed as follows: In the formula: For fluid diffusion residuals, For solid equilibrium residuals; The pressure is a dimensionless formation fluid pressure. For dimensionless time, for Dimensionless planar spatial coordinates of direction for Dimensionless planar spatial coordinates of direction This represents the dimensionless vertical displacement of a solid. , where is the Biot coefficient, characterizing the contribution ratio of pore pressure to effective stress; This represents the ratio of material constants.

5. The high-precision PINN multi-field coupling simulation method according to claim 1, characterized in that, In S5, the two-stage optimization strategy is as follows: First, a global search is performed using the Adam optimizer to initially determine the parameter range; then, the L-BFGS quasi-Newton algorithm is switched to perform local high-precision fitting, and the penetration rate is monitored. Elastic modulus By observing the evolutionary trajectory of the iteration steps, the true values ​​of permeability and elastic modulus are finally obtained.

6. A high-precision PINN multi-field coupling simulation system, characterized in that, Performing a high-precision PINN multi-field coupling simulation method as described in any one of claims 1-5 includes: The input module is used to acquire spatiotemporal coordinate data and initial guesses of the parameters to be inverted in the formation; The preprocessing module is used to perform dimensionless scaling preprocessing on the spatiotemporal coordinate data and the physical output quantities to be solved. The network configuration module is used to build a fully connected deep neural network and force the underlying operators and tensor data types to be declared as double-precision floating-point numbers. The residual construction module is used to construct the residual terms corresponding to the thermo-fluid-mechanical coupling law based on the porous elastic control equation using automatic differentiation technology. The sampling module is used to perform non-uniform encrypted sampling of the corner points and physical boundaries of the computational domain based on geometric feature perception. The loss function construction module is used to construct a multi-objective integrated loss function that includes physical constraints, observation data, and boundary conditions using an asymmetric weight configuration method. The parameter optimization module is used to set the parameters of the formation to be identified as learnable parameters of the network, and adopts a two-stage optimization strategy to collaboratively and iteratively update the network weights and the parameters to be identified. The output module is used to output the multiphysics simulation results and the identified formation permeability parameters after the iteration meets the preset convergence conditions.

Images