3D Holographic High-Performance Simulation Method and System Based on Multiscale Modeling Theory

Abstract

The present invention discloses a 3D holographic high-performance simulation method and system based on the multi-scale modeling theory. The method includes: realizing the adaptive meshing according to the change of wavelength through prior physical information to capture the details of the microscopic tissues, complex organs and multi-systems inside the human brain. At the same time, the wave velocity is incorporated into the adaptive mesh generation, overcoming the limitation of water-tight geometry in traditional CAD software. In order to avoid the overwhelming computational cost in high-frequency wave simulation, an efficient mass lumping method is adopted under the high-performance computing framework, and the numerical simulation of the two-dimensional and three-dimensional human brain is realized by combining triangular and tetrahedral meshes. The numerical method of the present invention has strong scalability and is suitable for the fast and accurate simulation of the full-brain full waveform. Through multi-scale whole-brain simulation, it can be used to study the influence of the skull on the propagation of high-frequency ultrasound in the brain, providing a new direction for the research and clinical application of ultrasound in the field of neurology.

Description

Technical Field

[0001] The present invention belongs to the field of biomedical simulation, and relates to a 3D holographic high-performance simulation method and system based on multi-scale modeling theory, which can be used to study the influence of the skull on the propagation of high-frequency ultrasound in the brain, and provides a new direction for the research and clinical application of ultrasound in the field of neurology medicine. Background Art

[0002] High-frequency ultrasound imaging plays a crucial role in brain medical research. Understanding the behavior of ultrasound in the brain at different scales from the system to the tissue is crucial for advancing high-resolution imaging. However, the complex internal structure of the brain poses significant challenges to precise geometric modeling, and due to the need to consider numerous wavelengths in space, the numerical simulation of high-frequency waves incurs a large amount of computational cost. Therefore, more efficient and accurate computational models and methods need to be developed. In the prior art, the Chinese invention patent "A Three-dimensional Modeling Method for Tomographic Image Anatomy of the Human Brain Striatum and Hippocampal Structure" (Publication No.: CN109544682A) discloses a three-dimensional modeling method for tomographic image anatomy of the human brain striatum and hippocampal structure. This method uses the three-dimensional fast spoiled gradient echo sequence of the MRI system to obtain high-resolution three-dimensional structure images of the whole brain in the horizontal plane, and then reconstructs the obtained images in the coronal plane and sagittal plane on a workstation, adjusts to the best contrast, observes simultaneously in the horizontal, coronal, and sagittal planes, and uses stereological methods for quantitative analysis to measure the relevant anatomical structures of the striatum and hippocampus; then uses SPSS13.0 statistical software to perform respective relevant statistical analyses on the data of the human brain striatum and hippocampal structure measured by MRI, and respectively compares and analyzes the data of the normal human brain striatum and hippocampal structure with the data under disease states; finally, imports the tomographic images into Mimics software to perform three-dimensional modeling on the data obtained by MRI. However, this method is limited by commercial software and cannot achieve the rapid generation of large-scale three-dimensional models. At the same time, the generated models are affected by the quality of MRI data, and there are problems such as insufficient local accuracy and modeling deviations caused by artifacts. Moreover, the generated CAD models are not suitable for finite element simulation.

[0003] To pursue high-resolution ultrasonic imaging and modeling of the human brain, on the one hand, it is necessary to generate high-quality mesh models, and on the other hand, advanced numerical calculation methods combined with high-performance computers are required to implement large-scale parallel computing. In mesh generation, it is simple and feasible to use a unified differential mesh for modeling (Guasch L, Calderón Agudo O, Tang M X, et al. Full-waveform inversion imaging of the human brain[J]. NPJ digital medicine, 2020, 3(1): 28.). However, due to the use of a unified mesh, this method faces accuracy limitations. In contrast, unstructured meshes provide flexibility for various material properties, facilitating simpler parameterization. For large-scale full-wave simulations of high-frequency waves, the finite element method based on high-performance computing is effective. However, when implementing an explicit time integration strategy, classical finite elements need to solve a linear system at each time step. The process of solving at each time step greatly reduces the computational efficiency.

[0004] Based on this, the present invention proposes a 3D full-brain holographic high-performance precise simulation technology based on the multi-scale modeling theory. It applies a mesh mapping method based on prior physical information to achieve spatial discretization of complex internal models of the human brain. At the same time, information such as frequency, solution time step, and mesh size range is incorporated into the prior information to constrain the generation of high-quality meshes. During the finite element numerical calculation process, the mass lumping technique is implemented to diagonalize the mass matrix to reduce the computational cost. At the same time, a two-level parallel strategy based on MPI communication is used for large-scale parallel computing to achieve precise modeling of the human brain model and full-wave numerical simulation of high-frequency ultrasonic waves. Since this method is an advanced numerical algorithm scheme based on a high-performance computing framework, it can be used for the full-wave modeling problem of any human brain model. Summary of the Invention

[0005] Aiming at the defects in the prior art, the purpose of the present invention is to provide a 3D full-brain holographic high-performance precise simulation method and system based on the multi-scale modeling theory.

[0006] The technical solution adopted by the present invention is as follows:

[0007] A 3D full-brain holographic high-performance precise simulation method based on the multi-scale modeling theory, including: based on the prior physical information of the human brain velocity field, combining the minimum wavelength change adjustment to achieve adaptive mesh refinement, using the generated mesh as a computational model, and solving high-frequency ultrasonic waves according to the residual form of the acoustic wave equation; the solution is to implement a full-wave numerical simulation of high-frequency ultrasonic waves using a high-order finite element algorithm based on mass-lumped elements under a two-level integrated parallel high-performance computing framework.

[0008] In the above technical solution, further, it specifically includes the following steps:

[0009] Step 1: According to the needs of simulation, select the prior physical information corresponding to the human brain velocity field, adjust the size of the edges of the constraint unit in combination with the minimum wavelength change, generate an adaptive mesh, and use it as a computational model.

[0010] Step 2: Add a perfectly matched layer to the model area to absorb the reflected waves from the boundary position.

[0011] Step 3: Implement the mass lumping technique in the solver, add bubble functions in the grid cells, and at the same time add high-order basis functions on the grid cells to diagonalize the mass matrix to reduce the computational time for inverting the mass matrix in each iteration, and obtain a mass lumped high-order finite element solver.

[0012] Step 4: Place the generated grid computational model and the established mass lumped high-order finite element solver under the high-performance computing framework, and realize the simulation by numerically solving the high-frequency ultrasonic waves through two-level integrated parallel computing.

[0013] Further, the prior physical information includes: the center frequency of the excitation source, the number of grid points per wavelength, the minimum time step, and the model velocity parameter.

[0014] Further, the generation of the adaptive mesh is specifically as follows:

[0015] The goal of optimizing the mesh is to maximize the numerical accuracy by adapting the mesh to the local changes in the physical field properties while minimizing the computational cost.

[0016] The mesh size is constrained by the CFL condition, that is:

[0017]

[0018] where h(x) is the diameter of the inscribed sphere or circle related to the spatial position element, dim is the spatial dimension, which is 2 or 3, is the simulated time step, is the local longitudinal sound wave velocity, is the Courant number which is a constant; thus, the minimum mesh size can be determined.

[0019] At the same time, the edge length of the model grid cell is related to the center frequency of the source wavelet and the longitudinal wave velocity of the sound wave:

[0020]

[0021] In the formula, G is the number of grid points per wavelength, is the sound wave wavelength, For each unit the degree of freedom, is the center frequency; at the same time, the unit edge length is related to the per-wavelength unit parameter C, , and the parameters C and G satisfy:

[0022]

[0023] wherein is a constant coefficient and is a function of the spatial polynomial k; thus, the grid size can be adjusted according to the minimum wavelength change.

[0024] Furthermore, the thickness of the perfectly matched layer set is greater than one quarter of the simulated minimum wavelength.

[0025] Furthermore, the solver uses an explicit time marching method to solve the second-order acoustic wave equation, and uses mass lumping to diagonalize the mass matrix to reduce the computational amount of inverse mass matrix calculation in each iteration. The diagonalization process includes: adding additional bubble functions in high-order grid elements, adding high-order basis functions on the grid elements, and using Lagrangian basis functions and Newton-Cauchy integration rules to achieve the diagonalization of the mass matrix.

[0026] Furthermore, the two-level integrated parallel computing realizes two-level integration parallelism of the excitation source and the computational model in time and space, specifically:

[0027] In the integrated parallel computing, it is assumed that n×n processes are called to compute the model, and each process communicates with each other through MPI. The communicator is split into many sub-communicators, and each sub-communicator is referred to as an integrated member. All sub-communicators are divided into two groups, one group is the integrated sub-communicators, and the other group is the source sub-communicators. For the integrated sub-communicators, the computational grid is divided into multiple regions, and each integrated member corresponds to a region, and finite element calculations of spatial parallelism are performed among the integrated members; for the source sub-communicators, each integrated member is equipped with an excitation source, and MPI communication is used among the integrated members to ensure the synchronization of time evolution during the computational execution process.

[0028] A computer-readable storage medium, on which a computer program is stored, and when the program is executed by a processor, it implements the 3D whole-brain holographic high-performance precise simulation method based on the multi-scale modeling theory as described in any one of the above.

[0029] A 3D whole-brain holographic high-performance precise simulation system based on the multi-scale modeling theory, comprising:

[0030] One or more processors;

[0031] A memory for storing one or more programs;

[0032] When the one or more programs are executed by the one or more processors, the one or more processors implement the 3D full-brain holographic high-performance precise simulation method based on the multi-scale modeling theory as described in any one of the above.

[0033] Compared with the prior art, the significant advantages of the present invention are as follows:

[0034] (1) An adaptive grid method that transcends the challenges of traditional CAD is proposed in the present invention: by adjusting the grid size according to the wavelength change, a high-quality hierarchical grid is generated using a grid generator, where the grid scale can be adjusted according to the change of the physical field characteristics of the internal tissues of the brain. High-fidelity, multi-scale human brain modeling is achieved. (2) An advanced numerical simulation method based on a high-performance computer is provided in the present invention: for large-scale simulation of high-frequency ultrasonic waves, with high-order accuracy and low computational cost, provided by an advanced mass lumping element, ensuring efficient and accurate simulation of wave propagation. (3) High accuracy: The numerical stability of high-order elements is ensured by applying an additional bubble function to the grid cells. At the same time, the accuracy of the solution is improved by using high-order basis functions. (4) Scalability: The use of mass lumping elements with a diagonal mass matrix eliminates the need to solve a linear system, thereby minimizing MPI communication and significantly improving parallel scalability. Brief Description of the Drawings

[0035] Figure 1 It is a schematic flowchart of the 3D full-brain holographic high-performance precise simulation method based on the multi-scale modeling theory of the present invention.

[0036] Figure 2 It is the parallel architecture under high-performance computing of the present invention: the construction of a space and source two-level integrated parallel processor.

[0037] Figure 3 It is the first application example of the present invention: a double-layer uniform small ball model.

[0038] Figure 4 It is the simulation result of the uniform small ball model: the pressure time series generated by the high-order mass lumping (k = 5) method and COMSOL at the receiver.

[0039] Figure 5 It is the signal difference between different solvers and the reference solution. The L2 norm errors of different degrees of high-order mass lumping methods (k = 3, 4, 5) and the reference solution COMSOL are 2.30%, 0.45%, and 0.15% respectively. The L2 norm error of COMSOL k = 5 and the reference solution is 0.18%.

[0040] Figure 6 It is the second application example of the present invention: a human brain adaptive grid distribution map based on prior physical information mapping. The grid contains a 1-inch thick PML.

[0041] Figure 7 It is the two - dimensional decomposition of the human brain region under a two - level integration framework.

[0042] Figure 8 It is the wave - field snapshots at different time steps with the skull.

[0043] Figure 9 It is the wave - field snapshots at different time steps without the skull.

[0044] Figure 10 It is the 50 - channel time - domain waveforms received with and without the skull (a. without skull; b. with skull).

[0045] Figure 11 It is the third application example of the present invention: three - dimensional non - invasive ultrasound experimental setup for the human brain.

[0046] Figure 12 In three - dimensional space, the brain adaptive grid distribution map based on prior physical information only shows local tissues or organs.

[0047] Figure 13 It is the three - dimensional decomposition of the human brain region under a two - level integration framework.

[0048] Figure 14 It is the wave - field snapshots at different time steps of the three - dimensional non - invasive ultrasound experiment for the human brain (sagittal plane).

[0049] Figure 15 It is the wave - field snapshots at different time steps of the three - dimensional non - invasive ultrasound experiment for the human brain (coronal plane).

[0050] Figure 16 It is the wave - field snapshots at different time steps of the three - dimensional non - invasive ultrasound experiment for the human brain (transverse plane). Detailed implementation manners

[0051] The technical solution of the present invention will be further described in detail below in conjunction with the accompanying drawings and specific examples. The following examples will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several changes and improvements can still be made. These all fall within the protection scope of the present invention.

[0052] The present invention proposes a 3D whole-brain holographic high-performance precise simulation technology based on multi-scale modeling theory, which is used to construct a complex multi-scale structure model of the human brain and simultaneously perform full-wave numerical simulation of high-frequency ultrasonic waves. The brain is filled with multi-scale structures, and huge grids are required to capture tiny tissues and complex organs. Aiming at the problem of modeling the complex internal structure of the human brain, the present invention uses prior physical information such as frequency, wavelength, time step, and speed to constrain the edges of the grid, and the generated grid cells adapt to the change of the shortest wavelength. Such grid cells with the shortest wavelength constraint can easily capture the tiny structures of the human brain in space. Aiming at the large amount of computational resource consumption caused by inverting the mass matrix at each time step in the explicit time evolution strategy, the present invention uses the mass lumping method to diagonalize the mass matrix, reducing the computational amount of inverting the mass matrix in each iteration. Aiming at the challenge of large-scale problem calculation, the present invention proposes a combination of a high-order finite element method based on the mass lumping method and a high-performance computing framework based on MPI two-level integrated parallelism, effectively utilizing computational resources and improving the computational efficiency of a wide range of problems.

[0053] According to a specific embodiment of the present invention, a 3D whole-brain holographic high-performance precise simulation method based on multi-scale modeling theory includes: based on the prior physical information of the human brain velocity field, adjusting and realizing adaptive grid meshing in combination with the change of the shortest wavelength, using the generated grid as a calculation model, and solving high-frequency ultrasonic waves according to the residual form of the acoustic wave equation; the solution is to implement full-wave numerical simulation of high-frequency ultrasonic waves using a high-order finite element algorithm based on mass lumped elements under a two-level integrated parallel high-performance computing framework.

[0054] Such as Figure 1 , the method specifically includes the following steps:

[0055] The first step: According to the needs of the simulation, select the corresponding prior physical information, including the center frequency of the excitation source, the number of grid points included in each wavelength, the minimum time step, the model speed parameter, etc., to constrain the size of the unit edges, and adjust the grid size according to the change of the shortest wavelength to generate an adaptive grid; this method can effectively capture the tiny tissues and complex organs inside the human brain, construct a multi-scale whole-brain model, and reduce the degrees of freedom required for calculation; this step can adaptively obtain two-dimensional triangular grids or three-dimensional tetrahedral grids according to the numerical simulation requirements of two-dimensional or three-dimensional human brains.

[0056] The second step: Add a perfectly matched layer to the model area to absorb the reflected waves from the boundary position; the thickness of the perfectly matched layer can be selected according to the simulated wavelength and the geometric size of the model, and preferably, the thickness of the set perfectly matched layer is greater than one-fourth of the simulated minimum wavelength.

[0057] Step 3: Implement the mass lumping technique in the solver to simplify the redundancy in the linear system caused by spatio-temporal discretization in wave problems, including adding additional bubble functions in grid cells to ensure the numerical stability of using high-order elements, and adding high-order basis functions on grid cells to reduce the computational time for inverting the mass matrix in each iteration while ensuring the solution accuracy; use Lagrangian basis functions and Newton-Cauchy integration rules to achieve the diagonalization of the mass matrix;

[0058] Step 4: Place the generated grid model and the established mass-lumped high-order finite element solver in a high-performance computing framework, and numerically solve the high-frequency ultrasonic waves through two-level integrated parallel computing, that is: Since large-scale parallel computing needs to be achieved, the excitation source and the model are integrated in parallel in time and space at two levels. By splitting the communicator into many sub-communicators and dividing them into two groups, each sub-communicator is regarded as an integration member. Each integration member in one group uses spatial parallelism in the sub-communicator in the usual way, and each integration member in the other group of sub-communicators is equipped with an excitation source, allowing communication between integration members. Achieve the full-wave simulation of the propagation of high-frequency ultrasonic waves inside the human brain.

[0059] The present invention can be used to reveal the influence of the skull on the propagation of high-frequency ultrasound in the brain, thereby promoting the research and clinical application of ultrasound in the field of neurology. The following details several key points in the solution of the present invention:

[0060] 1. The acoustic wave equation with the PML strategy applied

[0061] Ultrasonic waves are high-frequency acoustic waves that follow the wave equation of acoustic waves when propagating in space. Numerical simulation needs to simulate the space as an infinite region while reducing the contamination of the received signal due to boundary reflections. In the present invention, a perfectly matched layer is added as an absorbing boundary condition around the physical calculation region. Considering the propagation of the acoustic wave equation in the second-order form in a three-dimensional physical domain, the physical domain contains a wave absorption region wrapped by an external boundary PML layer and the model solution domain . The acoustic wave equation in the free radiation region can be expressed in the Cartesian coordinate system as:

[0062] (1)

[0063] In Equation (1), v is a free parameter related to the wave speed in space, called the longitudinal wave speed. The acoustic wave equation is approximately solved using the continuous Galerkin method, and the remaining quantity form is obtained through the introduction of auxiliary variables and Fourier transform.

[0064] (2)

[0065] In Equation (2), and respectively represent the introduced auxiliary vectors and auxiliary variables. Their expressions are

[0066] (3)

[0067] (4)

[0068] Suppose homogeneous Dirichlet boundary conditions are imposed on and , and their test functions are respectively , for no boundary condition is imposed. The first term in Equation (2) introduces an external excitation source f. In the computational domain, for time , at the moment, the coupled form of the modified second-order acoustic wave PML equations in the domain is obtained using the residual form as The coupled form of the modified second-order acoustic wave PML equations in the domain is

[0069] (5)

[0070] (6)

[0071] (7)

[0072] (8)

[0073] (9)

[0074] (10)

[0075] In the formula, represents the outer normal vector on the boundary , is the source vector, and the excitation source is a time-varying Ricker wavelet. is a Lipschitz boundary, is the unknown field. is the coefficient term related to the damping function .

[0076] 2. High-order mass lumping strategy

[0077] The explicit time marching method is used to solve the second-order acoustic wave equation. In order to avoid the large amount of computational resource consumption caused by the inverse solution of the mass matrix M at each time step in the explicit time marching strategy, mass lumping is adopted to diagonalize the mass matrix and reduce the computational amount of the inverse solution of the mass matrix in each iteration.

[0078] Equation (3-8) in the three-dimensional open computational physical domain with the damping function , the weak form of the acoustic wave equation under the action of a vector sound source can be expressed as:

[0079] (9)

[0080] Let be a tetrahedral tessellation that satisfies the general properties of the finite element mesh, and h be the diameter of the smallest inscribed circle of each mesh element within . The finite element space U h is defined as the space composed of continuous functions.

[0081] (10)

[0082] is a Hilbert space. represents the Lagrangian approximation subspace associated with . The problem of the classical conforming finite element method is to continuously find that satisfies:

[0083] (11)

[0084] Using to represent the sets of degrees of freedom and the corresponding Lagrangian bases associated with respectively. For any function , define as the coefficient vector. For any element, equation (11) can be rewritten using the linear basis h of U as

[0085] (12)

[0086] where are the mass matrix and the stiffness matrix respectively, expressed as

[0087] (13)

[0088] In the formula is the continuous inner product in L2.

[0089] Mass lumping is usually achieved using the nodal basis functions and the inexact orthogonality rules of the mass matrix. When the integration points coincide with the nodes of the basis functions, a diagonal matrix is obtained. During the diagonalization process, by adding additional bubble functions in high-order mesh elements and adding high-order basis functions on the mesh elements, the diagonalization of the mass matrix is achieved using the Lagrangian basis functions and the Newton-Cotes integration rules.

[0090] Consider the polynomial finite element space in the following form

[0091] (14)

[0092] New finite element subspace Satisfy . In the formula, , represents the space of high - order face and interior bubble functions. When k = 2, the new finite element consists of 7 Lagrangian interpolation points. There are three equivalence class points, including three vertices {S1, S2, S3}, the mid - points of three edges {M1, M2, M3} and the centroid point G.

[0093] Construct a new finite element subspace is to use a nodal basis , denotes the set of all nodes in The element nodes can be obtained by mapping from the reference node . Thus, by mapping the nodal basis to the physical element, the nodal basis functions of the element can be obtained

[0094] (15)

[0095] In the formula is the weight of the element . The global inner product can be expressed as:

[0096] (16)

[0097] According to the properties of the nodal basis functions, , represents the Kronecker function. Therefore, the corresponding mass matrix M ij can be expressed as

[0098] (17)

[0099] Since there are three equivalence classes of points in the element , the same weights are assigned to the points of the same class . The quadrature formula can be written as follows

[0100] (18)

[0101] For high - order finite elements, the quadrature formula can be derived in the same way. Compared with the numerical methods of the classical finite element method, the finite element method with lumped mass shows superior stability, convergence and efficiency.

[0102] 3. Adaptive Mesh Technology Based on Prior Physical Information Mapping

[0103] The following introduces the adaptive mesh technology to construct the finite element meshes required in the numerical calculation process. In this work, the goal of optimizing the mesh is to maximize the numerical accuracy by making the mesh adapt to the local changes in the physical field properties while minimizing the computational cost. To ensure numerical stability, the mesh size must be constrained by the CFL (Courant - Friedrichs - Lewey) condition because the time step is affected by the smallest element. Therefore, the mesh generation program needs to ensure that the smallest element is as large as possible to avoid overly small simulation time steps. For acoustic waves, the CFL condition that limits the mesh size can be described as

[0104] (19)

[0105] where \(h(x)\) is the diameter of the inscribed sphere or circle related to the spatial position element, dim is the spatial dimension of the problem (2 or 3), is the simulation time step, is the local acoustic longitudinal wave velocity. When the Courant number is defined as a constant (e.g., ), for a given and , rearrange Equation (19) to find the possible minimum mesh size.

[0106] Minimizing the computational cost requires the mesh to adapt to the local wavelength changes and adjust the mesh size according to the minimum wavelength changes; it is assumed that all triangular elements constituting the model are equilateral, which is necessary for accurate simulation using the finite element method. At this time, the edge length of the triangular element can be related to the center frequency of the source wavelet and the longitudinal wave velocity of the acoustic wave.

[0107] (20)

[0108] where \(G\) is the number of meshes per wavelength, is the acoustic wave wavelength, is the degree of freedom of each element , is the center frequency. At the same time, the edge length of the element can be related to the parameter \(c\) of the element per wavelength. Among them , the parameters \(C\) and \(G\) are related to each other in the following way:

[0109] (21)

[0110] where is a constant coefficient, which is a function of the spatial polynomial k. The mass lumped element has a higher number of nodes per element, so the higher-order mass lumped element has a higher value of each polynomial degree than the classical finite element.

[0111] 4. Two-Level Integrated Parallel Strategy Based on High-Performance Computing Framework

[0112] In integrated parallelism, it is assumed that n×n processes are called to compute the model. Each process communicates with each other through MPI, as Figure 2 shown. The communicator is split into many sub-communicators, and each sub-communicator is referred to as an integration member. All sub-communicators are divided into two groups, one group is the integration sub-communicators, and the other group is the source sub-communicators. For the integration sub-communicators, in each integration member, it is allowed to specify the existing functions of the finite element problem solution that uses spatial parallelism across sub-communicators in the normal way. In the other group of source sub-communicators, communication between integration members is allowed. By dividing the size of the original communicator by the number of processes in each integration member, each integration member will have the same spatial parallelism as the number of integration members. Therefore, the total number of processes initiated by mpiexec must be equal to the product of the number of set members and the number of processes used by each set member. In the actual calculation process, each integration member in the source sub-communicator is equipped with an excitation source, and MPI communication is used between integration members to ensure the synchronization of time evolution during the calculation execution. The integration members in the integration sub-communicator use sub-MPI communication to divide the computational grid into multiple regions.

[0113] Three numerical examples are given below:

[0114] (1) This case is to verify the accuracy and computational efficiency of the high-order finite element method using mass lumped elements. First, a two-layer uniform small ball model is set up, as Figure 3 shown. The wave velocity of the upper medium of the model is 1500 m / s, the wave velocity of the lower medium is 3200 m / s, and the wave velocity of the small ball is 2450 m / s. The size of the model area is 12-inch×12-inch, the excitation source coordinate is (6, 1.8), and the receiver coordinate is (6, 2.5). To verify the numerical accuracy, the model calculation uses high-order mass lumped finite elements of different degrees. The grid is encrypted by 10 times, and the solution obtained by the commercial simulation software COMSOL based on the discontinuous Galerkin finite element numerical method is used as the reference solution. Compare the time-domain waveforms obtained at the receiver between different methods, as Figure 4 shown.

[0115] When the program runs using different finite element solvers, the relationship between the degrees of freedom (DoF) and resource consumption is shown in Table 1. When using the high-order mass lumped ML, the numerical solution is approximated to the reference solution. However, the computational time scale is reduced by about 3 times.

[0116] Solver Degree of Freedom (DoF) Elapsed Time Memory (MB) ML k = 3 169848 9266 398 ML k = 4 287219 14789 443 ML k = 5 560818 10425 515 COMSOL k = 5 546798 31065 498 Reference 5172360 66324 4852

[0117] Table 1

[0118] Figure 5 The signal differences and L2 norm errors between the high-order mass lumping methods of different orders and the reference solution are given. The L2 error is calculated by the following formula:

[0119] (22)

[0120] where and represent the sound pressure of the numerical solution and the reference solution, respectively. The results show that the relative numerical error between the high-order mass lumping method and the reference solution decreases with the increase of the order, and the maximum L2 norm error does not exceed 0.25%, indicating that the high-order mass lumping method has good numerical stability when using high-order elements. The numerical method of the present invention has a high order of 5 and strong scalability, and is suitable for fast and accurate simulation of the whole brain full waveform.

[0121] (2) This case considers the experimental situation of simulating a human brain model placed in water. The speed of sound waves in water is set to a constant 1500 m / s. The size of the model area is 12-inch × 12-inch, and the speed of sound waves in each tissue inside the human brain varies in the range of 1500 m / s - 3000 m / s. A time-varying Ricker wavelet with a central frequency of 5.9 MHz is used as the excitation source, and an outer layer of unit grid with a thickness of 1-inch acts as the PML layer to absorb the reflected waves from the truncated boundary. The grid cell size is a function of the P-wave velocity of the sound wave as a variable, the constant G is set to 5, and mass lumping elements with k = 3 and a time step parameter of 0.0001 s are used to constrain the edge length of the elements. The computational domain is about 3600 wavelengths along one direction relative to the maximum cut-off frequency.

[0122] Figure 6 Shows the two-dimensional adaptive grid of the human brain generated by the mapping of prior physical information. Such grid cells with minimum wavelength constraints can easily capture the tiny structures of the human brain in space, but it has large cross-scale characteristics. For example, after locally magnifying 256 times in the cerebellar region, the changes in the speed of sound waves between the cerebellar cortex, medulla, and brainstem lead to scale changes in the mapped grid space. In addition, the adaptive grid conforms to the complex geometric structures of various brain tissues, and compared with the uniform structured grid, it reduces the degrees of freedom of the computational space. This method meets the accuracy requirements of specific solutions and adapts to local changes in physical properties.

[0123] Figure 7It presents a regional decomposition diagram of the two-dimensional human brain under a two-level integrated parallel high-performance computing framework. The model uses 80 cores for parallel computing and includes 10 excitation sources. Each excitation source divides the space into 8 regions for distributed memory parallel computing during the time evolution process. This method significantly improves the computing efficiency.

[0124] Figure 8 、 Figure 9 It shows the propagation process of ultrasound between the internal tissues of the human brain at different time steps calculated using high-order mass lumped elements combined with a two-level integrated parallel framework. The influence of having a skull ( Figure 8 ) and not having a skull ( Figure 9 ) on transcranial ultrasound is studied. Figure 10 It gives the time-domain waveforms received by the receiving array with and without a skull.

[0125] (3) This case further considers a non-invasive ultrasound experiment on the human brain in three-dimensional space. It simulates the scenario of placing a multi-modal human brain model generated by mapping prior physical information in a water tank. The size of the model area is 12-inch×12-inch×8.75-inch. The excitation source is placed on the right side of the water tank, and a 400×400 hydrophone array for receiving ultrasound signals is placed on the left, as Figure 11 shown. The model area is discretized using mass lumped elements with a degree of 3 (k = 3). The grid model generated using the adaptive grid technology based on prior physical information is as Figure 12 shown. The number of grid cells in the discretized area is 16450340, resulting in 114731163 degrees of freedom.

[0126] Under the two-level integrated parallel framework, the software runs on a server cluster consisting of six computing nodes. The server cluster includes 4 computing nodes with 56-core Intel Gold 6348 cpus, 1 112-core computing node, and 1 24-core management node. 300 cores are actually used during the execution process for parallel computing across multiple nodes. The regional decomposition of the model during the running process is as Figure 13 shown.

[0127] Figure 14 、 Figure 15 、 Figure 16 It shows the wave field snapshots of high-frequency ultrasound in the human brain at different time steps. By providing wave field snapshots of the sagittal plane, coronal plane, and transverse plane, it shows the propagation of ultrasound between the three-dimensional human brain tissues at different time points. These wave field snapshots provide valuable insights into the propagation laws and behavioral characteristics of ultrasound in a complex brain environment.

[0128] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk memory, CD-ROM, optical memory, etc.) that contain computer-usable program code.

[0129] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to the embodiments of the present invention. It should be understood that each flow and / or block in the flowchart and / or block diagram, and the combination of flows and / or blocks in the flowchart and / or block diagram, can be realized by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, such that the instructions executed by the processor of the computer or other programmable data processing devices generate means for realizing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0130] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, such that the instructions stored in the computer-readable memory generate a manufactured article including instruction means that realizes the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0131] These computer program instructions can also be loaded onto a computer or other programmable data processing device, such that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process, and thus the instructions executed on the computer or other programmable device provide steps for realizing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0132] The specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention.

Claims

1. A 3D holographic high-performance simulation method based on multi-scale modeling theory, characterized in that, Comprising: Based on the prior physical information of the human brain velocity field, adaptive mesh refinement is achieved by combining the adjustment of the minimum wavelength change. The generated mesh is used as a computational model, and the high-frequency ultrasonic waves are solved according to the residual form of the acoustic wave equation. The solution is realized by using a high-order finite element algorithm based on mass-lumped elements under a two-level integrated parallel high-performance computing framework for the full-wave numerical simulation of high-frequency ultrasonic waves. Among them, generating the adaptive mesh specifically includes: The goal of optimizing the mesh is to maximize the numerical accuracy by making the mesh adapt to the local changes of the physical field properties while minimizing the computational cost. The mesh size is constrained by the CFL condition, that is: ， where h(x) is the diameter of the inscribed sphere or circle related to the spatial position element, dim is the spatial dimension, which is 2 or 3, is the simulation time step, is the local acoustic longitudinal wave velocity, is the Courant number; from this, the minimum grid size is determined; Meanwhile, the edge length of the model grid cell is related to the center frequency of the source wavelet and the longitudinal wave velocity of the acoustic wave: ， where G is the number of wavelength grids per unit, is the acoustic wavelength, is the degree of freedom per element ; and is the center frequency; meanwhile, the edge length of the element is related to the per-wavelength element parameter C, , and the parameters C and G satisfy: ， wherein is a constant coefficient and is a function of the spatial polynomial k; thus, the grid size can be adjusted according to the minimum wavelength variation.

2. The 3D holographic high-performance simulation method based on the multi-scale modeling theory according to claim 1, wherein Specifically, it includes the following steps: The first step: According to the needs of the simulation, select the prior physical information corresponding to the human brain velocity field, combine the adjustment of the minimum wavelength change to constrain the size of the unit edges of the mesh, generate an adaptive mesh, and use it as a computational model. The second step: Add a perfectly matched layer to the model area to absorb the reflected waves from the boundary position. The third step: Implement the mass-lumping technique in the solver, add bubble functions in the mesh elements, and at the same time add high-order basis functions on the mesh elements to diagonalize the mass matrix to reduce the computational time for inverting the mass matrix in each iteration, and obtain a mass-lumped high-order finite element solver. The fourth step: Place the generated mesh computational model and the established mass-lumped high-order finite element solver under a high-performance computing framework, and numerically solve the high-frequency ultrasonic waves through two-level integrated parallel computing to achieve the simulation.

3. The 3D holographic high-performance simulation method based on the multi-scale modeling theory according to claim 2, characterized in that, The prior physical information described above includes: the center frequency of the excitation source, the number of mesh points per wavelength, the minimum time step, and the model velocity parameter.

4. The 3D holographic high-performance simulation method based on the multi-scale modeling theory according to claim 2, characterized in that, The thickness of the set perfectly matched layer is greater than one-fourth of the simulated minimum wavelength.

5. The 3D holographic high-performance simulation method based on the multi-scale modeling theory according to claim 2, characterized in that, The solver uses an explicit time marching method to solve the second-order acoustic wave equation, and uses mass-lumping to diagonalize the mass matrix to reduce the computational amount of inverse solving of the mass matrix in each iteration. The diagonalization process includes: adding additional bubble functions in the high-order mesh elements and adding high-order basis functions on the mesh elements, and using Lagrangian basis functions and Newton-Cauchy integration rules to diagonalize the mass matrix.

6. The 3D holographic high-performance simulation method based on the multi-scale modeling theory according to claim 2, wherein The two-level integrated parallel computing is to achieve two-level integration and parallelism of the excitation source and the computational model in time and space. Specifically: In the integrated parallel computing, assume that n×n processes are called to calculate the model. Each process communicates with each other through MPI. The communicator is split into many sub-communicators, and each sub-communicator is referred to as an integration member. All sub-communicators are divided into two groups, one group is the integration sub-communicator, and the other group is the source sub-communicator. For the integration sub-communicator, the computational mesh is divided into multiple regions, and each integration member corresponds to a region, and the integration members perform finite element calculations with spatial parallelism; for the source sub-communicator, each integration member is equipped with an excitation source, and MPI communication is used between the integration members to ensure the synchronization of the time evolution during the calculation execution.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by a processor, it implements the 3D holographic high-performance simulation method based on the multi-scale modeling theory described in any one of claims 1-6.

8. A 3D holographic high-performance simulation system based on the multi-scale modeling theory, characterized in that, Comprising: One or more processors; A memory for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors implement the 3D holographic high-performance simulation method based on the multi-scale modeling theory according to any one of claims 1-6.

Images