CN101770542A - Method of cluster computer for simulating electromagnetic wave propagation - Google Patents

Abstract

The invention relates to a method of a cluster computer for simulating electromagnetic wave propagation. The method comprises the following steps of: dispersing an electromagnetic wave propagation space into three-dimensional grids, dividing a simulation area into blocks by using an overlapping area disassembly method and constructing model data files; reading the corresponding model data files by each computation node of the cluster computer respectively, simulating the electromagnetic wave propagation by using a parallel electromagnetic field time-domain pseudo-spectrum method and outputting a file recording the steady-state distribution of an electromagnetic field in a space; and obtaining the distribution of a far-field scattered field by utilizing a near-field extrapolation far-field program. Compared with the traditional method, the invention has the following advantages that: the method can be used for simulating the electromagnetic wave propagation of point-source excitation; the large-scale electromagnetic wave propagation can be efficiently simulated in parallel; the iteration frequency in the simulation is reduced; and the far-field result is more accurate.

Description

The method of cluster computer for simulating electromagnetic wave propagation

Technical field

The present invention relates to electromagnetic wave propagation, particularly a kind of method of cluster computer for simulating electromagnetic wave propagation.

Background technology

Utilize numerical method that electromagnetic wave propagation is accurately simulated, most important to the rule of understanding electromagnetic wave and various complicated medium effects with the construction cycle that shortens electromagnetic product.Existing analogy method mainly contains method of moment, the pseudo-spectral method of electromagnetic field Finite Difference-Time Domain separating method and electromagnetic field time domain etc.

In the mentioned in the above the whole bag of tricks, the pseudo-spectral method of electromagnetic field time domain in 1997 by Q.H.Liu propose (Microwave and optical technology letters, volume: 15, page or leaf: 158-165,1997), relatively be applicable to the simulation of large-scale electromagnetic wave propagation.It is 3D grid with the spatial division of electromagnetic wave propagation that this method is at first chosen suitable grid length, is the propagation of analog electrical magnetic wave at infinite space, the outermost n of 3D grid _AbsorbLayer is an electromagnetic wave absorbing layer.Electric field and magnetic field are positioned at the summit of each grid cell in the space, draw time t by time iteration then to the Maxwell equation group after electromagnetic field in spatial distributions.Obtain the distribution of far field scattered field at last by extrapolation far field, near field program.In simulation process, space derivative in the Maxwell equation group is calculated by the pseudo-spectral method of a kind of high-precision Fourier, according to Nyquist's theorem, grid length in space usually only need be slightly less than electromagnetic wave half of minimal wave length when propagating in each medium of space in this method, thereby can greatly reduce the number of variable in the analog computation.

Yet there is the defective of the following aspects in the pseudo-spectral method of existing electromagnetic field time domain when using:

(1) electric field and magnetic field overlap on the locus in this method, and the space derivative operator that is adopted under this building method can produce Gibbs phenomenon when acting on discontinuous function.Thereby the structure of existing method only can be simulated the electromagnetic wave propagation problem of extended source excitation well.For the electromagnetic wave propagation simulation of a source forcing, need to adopt the spatial spread source to be similar to usually and replace the some source forcing, this way can not definitely reflect the real simulation environment.

(2) in a lot of electromagnetic wave propagations simulations (for example light wave in biological tissue, propagate simulation), the number of variable is quite huge, and the memory size of single processor, and simulated time head and shoulders above also can exceed acceptable scope.In order to address this problem effectively, can utilize cluster computer to carry out parallel processing usually, whole zoning is decomposed each computing node of cluster computer, each computing node calculates the propagation condition of electromagnetic field in the region separately respectively then.Yet because used space derivative operator is an overall operator in the pseudo-spectral method of electromagnetic field time domain, it need use the electromagnetic field data in the whole zone, and this just exchanges mass data inevitably between each computing node of whole cluster computer.This fact make the pseudo-spectral method of existing electromagnetic field time domain can't be on the distributed memory cluster computer parallel processing efficiently with the simulation large-scale electromagnetic wave propagation.

(3) time diffusion of the pseudo-spectral method of electromagnetic field time domain calculates the difference scheme that adopts second order accuracy, when time step length is not can cause bigger numerical dispersion error in enough hour.This point is particularly serious for large-scale electromagnetic field simulation problem, because the phase error that the numerical dispersion error is brought can constantly accumulate along with the increase of electromagnetic wave propagation distance.Usually need need choose in calculating the time step much smaller than stability condition, this has increased the needed iteration step number of simulation undoubtedly.

(4) lattice point that is adopted in the pseudo-spectral method of electromagnetic field time domain is coarse, thereby bigger at the enterprising line number value of sealing surface integration time error to equivalent surface current in obtaining extrapolation far field, the near field program of far field electromagnetic field distribution.

Therefore,, must improve, need a kind of efficient feasible parallel scheme especially to solve the large-scale electromagnetic wave propagation problem to existing method in order to adapt to the simulation of different electromagnetic wave propagation problems.

Summary of the invention

Fundamental purpose of the present invention provides a kind of method of cluster computer for simulating electromagnetic wave propagation, makes it can simulate the large-scale electromagnetic wave propagation problem efficiently.The present invention also is intended to expand the scope of application of the pseudo-spectral method of electromagnetic field time domain, improves the far field precision of calculation results.

Technical scheme of the present invention is as follows:

A kind of method of cluster computer for simulating electromagnetic wave propagation is characterized in that, this method comprises the steps:

(1) the electromagnetic wave propagation problem of simulation is set up model and model of creation data file as required:

1. needing the electromagnetic wave propagation problem simulated for the electromagnetic wave of certain form incides behind the area of space Ω in the propagation in this zone, is certain single medium outside the regional Ω; Choose three mutually orthogonal directions arbitrarily With

Satisfy

With Have the rectangular parallelepiped A of the minimum of an inclusion region Ω, each limit of A is parallel to Or

Center with A is an initial point,

With The positive dirction that is respectively x axle, y axle and z axle is set up three-dimensional cartesian coordinate system;

2. choose grid cell length Δ x, Δ y and Δ z;

3. be 3D grid according to selected grid cell length with spatial division, be respectively N ' at the number of x, y and z direction grid cell _x, N ' _yAnd N ' _z, guarantee that rectangular parallelepiped A takes up space in 3D grid; Determine to meet following formula and meet and to make the fast fourier transformation algorithm minimum positive integer L of characteristic length efficiently _x, L _yAnd L _z:

L _x≥((N′ _x+2n _absorb+20)+2n _overlap(n _x-1))/n _x

L _y≥((N′ _y+2n _absorb+20)+2n _overlap(n _y-1))/n _y

L _z≥((N′ _z+2n _absorb+20)+2n _overlap(n _z-1))/n _z

N wherein _OverlapBe the thickness of overlapping region in the model, n _AbsorbBe the thickness of electromagnetic wave absorbing layer, n _x, n _yAnd n _zBe respectively the number of cluster computer three-dimensional topology structure at x, y and z direction calculating node; Calculate

N _x＝2(L _x-n _overlap)+(n _x-2)(L _x-2n _overlap)

N _y＝2(L _y-n _overlap)+(n _y-2)(L _y-2n _overlap)

N _z＝2(L _z-n _overlap)+(n _z-2)(L _z-2n _overlap)；

4. with existing 3D grid at x axle positive dirction and the x axle negative direction n that stretches out respectively _X1And n _X2Layer grid cell obtains a new 3D grid, and the number of x, y and z direction grid cell is respectively N _x, N ' _yAnd N ' _z, wherein

n _X2=N _x-N ' _x-n _X1, function ceil (s) gets the minimum positive integer more than or equal to independent variable s;

5. with the 4. the new 3D grid that obtains of step at y axle positive dirction and the y axle negative direction n that stretches out respectively _Y1And n _Y2Layer grid cell obtains a new 3D grid, and the number of x, y and z direction grid cell is respectively N _x, N _yAnd N ' _z, wherein

n _Y2=N _y-N ' _y-n _Y1

6. with the 5. the new 3D grid that obtains of step at z axle positive dirction and the z axle negative direction n that stretches out respectively _Z1And n _Z2Layer grid cell obtains final 3D grid, and the number of x, y and z direction grid cell is respectively N _x, N _yAnd N _z, wherein

n _Z2=N _z-N ' _z-n _Z1The occupied space of this 3D grid is the whole space of electromagnetic wave propagation simulation, outermost n _AbsorbLayer grid cell is electromagnetic wave absorbing layer;

7. the method according to the Yee grid disposes electric field and the position of magnetic field in each grid cell of 3D grid; Removal is positioned at all electromagnetic field components on positive dirction one side end face of 3D grid all directions, respectively corresponding three-dimensional matrice of each electromagnetic field component in the grid, and these matrixes are of similar shape: the length of x direction is N _x, the length of y direction is N _y, the length of z direction is N _z

8. the whole spatial division of utilizing the overlapping region decomposition method to simulate is piece, the number of the piece of x, y and z direction equals the number of this direction calculating node in the cluster computer three-dimensional topology structure respectively, the scope of each electromagnetic field component is determined by following method in each piece: the matrix that each electromagnetic field component is formed is divided into continuous submatrix respectively, the number of all directions submatrix equals the number of the piece of this direction respectively, and the number of each submatrix element on direction ξ is L _ξ-n _OverlapOr L _ξ-2n _Overlap, when this submatrix is the head of ξ direction or odd amount in addition to the round number matrix, be L _ξ-n _Overlap, all the other situations are L _ξ-2n _Overlap, wherein ξ is x, y or z; The zone of each submatrix of gained is extended, principle is if this submatrix has in adjacent submatrix and this submatrix corresponding electromagnetic field component direction and this border tangent on certain border, then with this submatrix at this boundary coordinate of border vertical direction n that stretches out therewith _OverlapEach submatrix of above-mentioned each electromagnetic field component is one by one corresponding to each piece, and the zone before extending is the Non-overlapping Domain in this piece, and the zone of Yan Shening was the overlapping region afterwards;

9. access time, step delta t was

Wherein T is the electromagnetic cycle, v _MaxThe maximum rate of in each medium of space, propagating for electromagnetic wave;

10. determine the time step n that iterations N and record electromagnetic field distribute ₁And n ₂

The specific inductive capacity and the magnetic permeability of each medium in the space are revised as ε respectively _h'=α ε _hAnd μ _h'=α μ _h, wherein

Be electromagnetic angular frequency, ε _hAnd μ _hBe respectively the specific inductive capacity and the magnetic permeability of h class medium reality, ε _h' and μ _h' be respectively the h class medium specific inductive capacity and the magnetic permeability that in simulation, use;

The model of creation data file: described overlapping region decomposition method is cut apart each piece that obtains, and sets up the model data file respectively, and called after " model_i _b_ j _b_ k _b", i wherein _b, j _bAnd k _bBe respectively the sequence number of this piece, following content write the pairing model data file of each piece in x, y and z direction: in this piece each electromagnetic field component in the coordinate range of all directions overlapping region and Non-overlapping Domain, this piece in the Non-overlapping Domain each amended specific inductive capacity of electromagnetic field component position medium and magnetic permeability, grid cell at the length of all directions, time step, iterations, the time step n in the simulation ₁And n ₂, and the parameter of incident field;

(2) cluster computer carries out data processing:

Read the model data file, with the pseudo-spectral method Parallel Simulation of parallel electromagnetic field time domain electromagnetic wave propagation, output record time step n ₁And n ₂The destination file that the time space electromagnetic field distributes;

1. the model data file of resulting correspondence in each computing node difference read step (1) of cluster computer obtains simulating required parameter, and coordinate is (i in the cluster computer three-dimensional topology structure _p, j _p, k _p) the model data file of computing node correspondence be " model_i _p_ j _p_ k _p"; The space of each computing node each electromagnetic field component of memory allocated in internal memory;

2. utilize the electromagnetic wave propagation in the pseudo-spectral method virtual space of parallel electromagnetic field time domain;

3. export destination file " result1 " and " result2 ", this destination file " result1 " and " result2 " have write down n respectively ₁Step and n ₂The numerical value of step each electromagnetic field component of time space;

(3) scattered field that utilizes extrapolation far field, near field program to calculate the far field distributes:

1. destination file " result1 " and " result2 " of output in the read step (2) obtain the complex field of each electromagnetic field component in the near field;

2. in 3D grid, the sealing surface of a virtual enclosing region Ω is set, by six tangent plane x=x ₁, x=x ₂, y=y ₁, y=y ₂, z=z ₁And z=z ₂Surround;

3. interpolation obtains each vertical with this tangent plane normal direction on above-mentioned each tangent plane magnetic-field component;

4. on each tangent plane, respectively each equivalent face electric current and each equivalent face magnetic current are made two-dimension fourier transform, each the equivalent face electric current on each tangent plane and each equivalent face magnetic current are expanded into the form of sum of series;

5. resolve integration and obtain the direction of extrapolating

The far field time required intermediate quantity

With

6. by With

Calculate direction

The far field scattered field

With

The pseudo-spectral method of described parallel electromagnetic field time domain is by time iteration simulation electromagnetic wave propagation, and electromagnetic field is finished following four steps successively by step current time iteration each computing node to the process of next time step in the space:

(1) utilize the time iteration formula of electric field to ask the electric field of Non-overlapping Domain in next this computing node of time step; Wherein the space derivative of each magnetic-field component that relates in the time iteration formula of electric field calculates with the pseudo-spectral method of the Fourier of this intranodal, promptly only operate at the magnetic field data of this computing node, comprise its overlapping region and Non-overlapping Domain, and earlier on the magnetic field of overlapping region, be multiplied by a smooth function and make it drop to 0 smoothly from the end to end that is connected with Non-overlapping Domain;

(2) receive the electric field data of the overlapping region that each neighborhood calculation node sends over, and send the electric field data that the overlapping region of self Non-overlapping Domain and this neighborhood calculation node partially overlaps to each neighborhood calculation node;

(3) utilize the time iteration formula in magnetic field to ask the magnetic field of Non-overlapping Domain in next this computing node of time step; Wherein the space derivative of each electric field component that relates in the time iteration formula in magnetic field calculates with the pseudo-spectral method of the Fourier of this intranodal, promptly only operate at the electric field data of this computing node, comprise its overlapping region and Non-overlapping Domain, and earlier on the electric field of overlapping region, be multiplied by a smooth function and make it drop to 0 smoothly from the end to end that is connected with Non-overlapping Domain;

(4) receive the magnetic field data of the overlapping region that each neighborhood calculation node sends over, and send the magnetic field data that the overlapping region of self Non-overlapping Domain and this neighborhood calculation node partially overlaps to each neighborhood calculation node.

Technique effect of the present invention:

(1) in the pseudo-spectral method of electromagnetic field time domain,, can not produce Gibbs phenomenon when acting on point source based on the space derivative operator of this building method with the position in electric field in the 3D grid and the magnetic field grid length that on all directions of space, staggers half.The electromagnetic wave propagation that building method after the improvement not only can be used for accurately simulating the electromagnetic wave propagation of extended source excitation but also can be used for accurate simulation points source forcing.

(2) adopt the overlapping region decomposition method that the pseudo-spectral method of electromagnetic field time domain is carried out parallel processing.Whole zoning is divided into piece and is assigned in each computing node of cluster computer and go, the data of each computing node are divided into overlapping region and Non-overlapping Domain.In the computation process of each time step, each node only calculates the electromagnetic field of Non-overlapping Domain separately, and the electromagnetic field of overlapping region is obtained by adjacent node by exchanges data after the calculating of each time step.The particularly important is, when calculating the electromagnetic field of Non-overlapping Domain, only need utilize data in each computing node internal memory to come the derivative of computer memory electromagnetic field.In order in this process, not cause error, the overlapping region part in the local data is multiplied by smooth function to form taper end.This method for parallel processing has been avoided exchanging mass data in each computing node of cluster computer, and exchanges data only limits to a spot of overlapping region between the neighborhood calculation node, thereby can simulate large-scale electromagnetic wave propagation efficiently.

(3) specific inductive capacity and the magnetic permeability of each medium of space that is adopted in the modification simulation.By in simulation, revising the specific inductive capacity and the magnetic permeability of each medium of space, the present invention has eliminated the bigger numerical dispersion error of being brought when choosing bigger time step, thereby increase adoptable time step in the simulation exponentially, reduce the required iteration step number of simulation, save computational resource.

(4) equivalent surface current and the integration method of equivalent face magnetic current on sealing surface in the extrapolation far field calculating of near field are improved.At first equivalent face electric current on the extrapolated side and equivalent face magnetic current are made two-dimension fourier transform acquisition fourier spectrum, thereby they can be expanded into the form of sum of series, partly can resolve integration for the fast-changing phase delay in space.This improvement has overcome the coarse shortcoming of lattice point in the pseudo-spectral method of electromagnetic field time domain, and far field result of calculation is more accurate.

In a word, advantage of the present invention is the electromagnetic wave propagation that can be used for the simulation points source forcing; Parallel Simulation large-scale electromagnetic wave propagation efficiently; Reduce number of iterations in the simulation; The far field result is more accurate.

Description of drawings

Fig. 1 is that space lattice of the present invention is divided synoptic diagram

Fig. 2 is electric field on each grid cell and a magnetic field position synoptic diagram in the space lattice of the present invention

Fig. 3 is that cluster computer of the present invention reads model data file, transaction module and exports the schematic flow sheet of destination file

Fig. 4 is the synoptic diagram of the smooth function in the pseudo-spectral method of the parallel electromagnetic field time domain of the present invention

Fig. 5 is the synoptic diagram that each computing node is asked the electromagnetic field space derivative in the pseudo-spectral method of the parallel electromagnetic field time domain of the present invention

Fig. 6 is the synoptic diagram of overlapping region of the present invention decomposition method

Fig. 7 is the synoptic diagram of the embodiment of the invention 1 midplane ripple incident biological tissue model

Fig. 8 is the schematic flow sheet of the embodiment of the invention 1

Fig. 9 is the analog result of the embodiment of the invention 1

Figure 10 is the schematic flow sheet of the embodiment of the invention 2

Figure 11 is the analog result of the embodiment of the invention 2

Embodiment

The invention will be further described below in conjunction with drawings and Examples, but should not limit protection scope of the present invention with this.

The method of cluster computer for simulating electromagnetic wave propagation among the present invention comprises three steps:

One, set up model and model of creation data file according to the electromagnetic wave propagation problem of simulation:

1. needing the electromagnetic wave propagation problem simulated for the electromagnetic wave of certain form incides behind the area of space Ω in the propagation in this zone, is certain single medium outside the regional Ω; Choose three mutually orthogonal directions arbitrarily

With

Satisfy

With

Have the rectangular parallelepiped A of the minimum of an inclusion region Ω, each limit of A is parallel to

Or

Center with A is an initial point,

With The positive dirction that is respectively x axle, y axle and z axle is set up three-dimensional cartesian coordinate system;

2. choose grid cell length Δ x, Δ y and Δ z; The system of selection of all directions grid cell length is consistent with the pseudo-spectral method of existing electromagnetic field time domain, must less than electromagnetic wave in each medium of space minimal wave length 1/2, on this basis must not be greater than the uneven smallest dimension of this direction electromagnetic property;

3. be 3D grid according to selected grid cell length with spatial division, be respectively N ' at the number of x, y and z direction grid cell _x, N ' _yAnd N ' _z, guarantee that rectangular parallelepiped A takes up space in 3D grid; Determine to meet following formula and meet and to make the fast fourier transformation algorithm minimum positive integer L of characteristic length efficiently _x, L _yAnd L _z:

L _x≥((N′ _x+2n _absorb+20)+2n _overlap(n _x-1))/n _x

L _y≥((N′ _y+2n _absorb+20)+2n _overlap(n _y-1))/n _y

L _z≥((N′ _z+2n _absorb+20)+2n _overlap(n _z-1))/n _z

N wherein _OverlapBe the thickness of overlapping region in the model, n _AbsorbBe the thickness of absorption layer, n _x, n _yAnd n _zBe respectively the number of cluster computer three-dimensional topology structure at x, y and z direction calculating node, can make fast fourier transform efficiently characteristic length l satisfy form l=2 ^m3 ⁿ5 ^d(m, n and d are nonnegative integer); Calculate

N _x＝2(L _x-n _overlap)+(n _x-2)(L _x-2n _overlap)

N _y＝2(L _y-n _overlap)+(n _y-2)(L _y-2n _overlap)

N _z＝2(L _z-n _overlap)+(n _z-2)(L _z-2n _overlap)；

4. with existing 3D grid at x axle positive dirction and the x axle negative direction n that stretches out respectively _X1And n _X2Layer grid cell obtains a new 3D grid, and the number of x, y and z direction grid cell is respectively N _x, N ' _yAnd N ' _z, wherein

n _X2=N _x-N ' _x-n _X1, function ceil (s) gets the minimum positive integer more than or equal to independent variable s;

5. with the 4. the new 3D grid that obtains of step at y axle positive dirction and the y axle negative direction n that stretches out respectively _Y1And n _Y2Layer grid cell obtains a new 3D grid, and the number of x, y and z direction grid cell is respectively N _x, N _yAnd N ' _z, wherein

n _Y2=N _y-N ' _y-n _Y1

6. with the 5. the new 3D grid that obtains of step at z axle positive dirction and the z axle negative direction n that stretches out respectively _Z1And n _Z2Layer grid cell obtains final 3D grid as shown in Figure 1, and the number of x, y and z direction grid cell is respectively N _x, N _yAnd N _z, wherein

n _Z2=N _z-N ' _z-n _Z1The occupied space of this 3D grid is the whole space of electromagnetic wave propagation simulation, outermost n _AbsorbLayer grid cell is electromagnetic wave absorbing layer;

7. according to the method configuration electric field and the position of magnetic field in 3D grid of Yee grid, as shown in Figure 2: electric field is positioned at the mid point of each seamed edge of grid cell, and direction of an electric field is along the positive dirction of place seamed edge; Magnetic field is positioned at the center of each grid cell face, and magnetic direction is along the positive dirction of the vertical line of place face; Removal is positioned at all electromagnetic field components on positive dirction one side end face of 3D grid all directions, respectively corresponding three-dimensional matrice of each electromagnetic field component in the grid, and these matrixes are of similar shape, and the length of x direction is N _x, the length of y direction is N _y, the length of z direction is N _z

8. the whole spatial division of utilizing the overlapping region decomposition method to simulate is piece, the number of the piece of x, y and z direction equals the number of this direction calculating node in the cluster computer three-dimensional topology structure respectively, the scope of each electromagnetic field component is determined by following method in each piece: the matrix that each electromagnetic field component is formed is divided into continuous submatrix respectively, the number of all directions submatrix equals the number of piece on this direction respectively, and the number of each submatrix element on direction ξ is L _ξ-n _OverlapOr L _ξ-2n _Overlap, when this submatrix is the head of ξ direction or odd amount in addition to the round number matrix, be the value in front, all the other situations are a next value, wherein ξ is x, y or z; The zone of each submatrix of gained is extended, method is if this submatrix has in adjacent submatrix and this submatrix corresponding electromagnetic field component direction and this border tangent on certain border, then with this submatrix at this boundary coordinate of border vertical direction n that stretches out therewith _OverlapEach submatrix of above-mentioned each electromagnetic field component that obtains is one by one corresponding to each piece, and the zone before extending is the Non-overlapping Domain in this piece, and the zone of Yan Shening was the overlapping region afterwards; Can provide the coordinate range of each electromagnetic field component in each piece by said method, with E _xFor example provides expression formula: because E _xAlong the x direction, therefore E on the x direction _xDo not have overlappingly, and between y and z direction and neighborhood calculation node, exist overlapping.E in i piece on the x direction _xData interval [a _i, b _i] can determine by following formula:

$[a_i,b_i]=\left\{\begin{array}{cc}[1,L_x-n_{\mathrm{overlap}}]& i=1\\ [b_{i-1}+1,a_i+L_x-2_{\mathrm{noverlap}}-1]& i=2,...,n_x-1\\ [b_{i-1}+1,N_x]& i=n_x\end{array}\right.$

E in the j piece on the y direction _xData interval [c _j, d _j] determine by following formula

$[c_j,d_j]=\left\{\begin{array}{cc}[1,L_y]& j=1\\ [d_{j-1}-2n_{\mathrm{overlap}}+1,c_j+L_y-1]& j=2,...,n_y-1\\ [d_{j-1}-2n_{\mathrm{overlap}}+1,N_y]& j=n_y\end{array}\right.$

E in the k piece on the z direction _xData interval [e _k, f _k] determine by following formula

$[e_k,f_k]=\left\{\begin{array}{cc}[1,L_z]& k=1\\ [f_{k-1}-2n_{\mathrm{overlap}}+1,e_k+L_z-1]& k=2,...,n_z-1\\ [f_{k-1}-2n_{\mathrm{overlap}}+1,N_z]& k=n_z\end{array}\right.$

The coordinate range of all the other each electromagnetic field components similarly;

9. access time, step delta t was

Wherein

T is the electromagnetic cycle, v _MaxThe maximum rate of in each medium of space, propagating for electromagnetic wave;

10. determine the iterations in the simulation

T wherein _SimulationBe the time of electromagnetic wave propagation in the simulation, also promptly in the space electromagnetic field distribute and reach the needed time of stable state, be chosen for the twice that electromagnetic wave comes and goes whole virtual space required time;

Time step n with electromagnetic field distribution in the record space ₁And n ₂Be set to n respectively ₁=N-m and n ₂=N;

The specific inductive capacity and the magnetic permeability of each medium in the space are revised as ε respectively _h'=α ε _hAnd μ _h'=α μ _h, wherein

Be electromagnetic angular frequency, ε _hAnd μ _hBe respectively the specific inductive capacity and the magnetic permeability of h class medium reality, ε _h' and μ _h' be respectively the h class medium specific inductive capacity and the magnetic permeability that in simulation, use;

The model of creation data file: cut apart each fritter that obtains for the overlapping region decomposition method, set up a model data file, and called after " model_i _b_ j _b_ k _b", i wherein _b, j _bAnd k _bBe respectively the sequence number of this piece in x, y and z direction; Following content is write the pairing model data file of each piece: in this piece each electromagnetic field component in the coordinate range of all directions overlapping region and Non-overlapping Domain, this piece in the Non-overlapping Domain each amended specific inductive capacity of electromagnetic field component position medium and magnetic permeability, grid cell at the length of all directions, time step, iterations, the time step n in the simulation ₁And n ₂And the parameter of incident wave; The parameter of incident wave includes all required parameters of each component value of electromagnetic field of definite incident wave any moment arbitrfary points in the space such as ejected wave shape, direction and frequency;

Two, cluster computer carries out data processing:

Read the model data file, with the pseudo-spectral method Parallel Simulation of parallel electromagnetic field time domain electromagnetic wave propagation, output record time step n ₁And n ₂The destination file that the time space electromagnetic field distributes;

Fig. 3 has showed the process flow diagram of this step:

(1) message call passing interface program starts MPI (Message Passing Interface) process;

Start p process (p=n _x* n _y* n _z), correspond respectively to p computing node in the cluster computer; This p process is mapped to the three-dimensional topology structure, and the number of x, y and z direction is respectively n _x, n _yAnd n _z

(2) each computing node reads pairing model file;

Each computing node reads corresponding model file respectively and (supposes that the coordinate of certain computing node in the three-dimensional topology structure is (i _p, j _p, k _p), then its pairing model file is " model_i _p_ j _p_ k _p"), obtain following parameter: length Δ x, Δ y and the Δ z of grid cell all directions, the time step n that electromagnetic field distributes in the time step Δ t in the simulation, total iterations N, the record space ₁And n ₂, each electromagnetic field component Non-overlapping Domain and overlapping region coordinate range, each electromagnetic field component in the parameter of the electromagnetic parameter and the incident wave of Non-overlapping Domain position; The memory headroom of each electromagnetic field component of memory allocated;

(3) begin simulation, time step n=1 is set;

(4) each computing node utilizes the time iteration formula of electric field to calculate in the Non-overlapping Domain next step electric field;

The time iteration formula of electric field is divided into two zones:

Non-absorption layer,

$\left\{\begin{array}{c}{E}_x{|}_{i+1/2,j,k}^{n+1}=E_x{|}_{i+1/2,j,k}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ}}^{\′}{|}_{i+1/2,j,k}}[(\frac{\∂H_z}{\∂y}-\frac{\∂H_y}{\∂z}){|}_{i+1/2,j,k}^{n+1/2}-({\mathrm{\ϵ}}^{\′}{|}_{i+1/2,j,k}-{\mathrm{\ϵ}}_0^{\′})\·\left(\frac{\∂E_{x,\mathrm{inc}}}{\∂t}\right){|}_{i+1/2,j,k}^{n+1/2}]\\ E_y{|}_{i,j+1/2,k}^{n+1}=E_y{|}_{i,j+1/2,k}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ}}^{\′}{|}_{i,j+1/2,k}}[(\frac{\∂H_x}{\∂z}-\frac{\∂H_z}{\∂x}){|}_{i,j+1/2,k}^{n+1/2}-({\mathrm{\ϵ}}^{\′}{|}_{i,j+1/2,k}-{\mathrm{\ϵ}}_0^{\′})\·\left(\frac{\∂E_{y,\mathrm{inc}}}{\∂t}\right){|}_{i,j+1/2,k}^{n+1/2}]\\ E_z{|}_{i,j,k+1/2}^{n+1}=E_z{|}_{i,j,k+1/2}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ}}^{\′}{|}_{i,j,k+1/2}}[(\frac{\∂H_y}{\∂x}-\frac{\∂H_x}{\∂y}){|}_{i,j,k+1/2}^{n+1/2}-({\mathrm{\ϵ}}^{\′}{|}_{i,j,k+1/2}-{\mathrm{\ϵ}}_0^{\′})\·\left(\frac{\∂E_{z,\mathrm{inc}}}{\∂t}\right){|}_{i,j,k+1/2}^{n+1/2}]\end{array}\right.$

Absorption layer,

$\left\{\begin{array}{c}{P}_x{|}_{i+1/2,j,k}^{n+1}=\frac{2-\mathrm{\Δt}\·b_y{|}_j}{2+\mathrm{\Δt}\·b_y{|}_j}{P}_x{|}_{i+1/2,j,k}^n+\frac{2\mathrm{\Δt}}{(2+\mathrm{\Δt}\·b_y{|}_j){\mathrm{\ϵ}}_0^{\′}}(\frac{\∂H_z}{\∂y}-\frac{\∂H_y}{\∂z}){|}_{i+1/2,j,k}^{n+1/2}\\ P_y{|}_{i,j+1/2,k}^{n+1}=\frac{2-\mathrm{\Δt}\·b_z{|}_k}{2+\mathrm{\Δt}\·b_z{|}_k}{P}_y{|}_{i,j+1/2,k}^n+\frac{2\mathrm{\Δt}}{(2+\mathrm{\Δt}\·b_z{|}_k){\mathrm{\ϵ}}_0^{\′}}(\frac{\∂H_x}{\∂z}-\frac{\∂H_z}{\∂x}){|}_{i,j+1/2,k}^{n+1/2}\\ P_z{|}_{i,j,k+1/2}^{n+1}=\frac{2-\mathrm{\Δt}\·b_x{|}_i}{2+\mathrm{\Δt}\·b_x{|}_i}{P}_z{|}_{i,j,k+1/2}^n+\frac{2\mathrm{\Δt}}{(2+\mathrm{\Δt}\·b_x{|}_i){\mathrm{\ϵ}}_0^{\′}}(\frac{\∂H_y}{\∂x}-\frac{\∂H_x}{\∂y}){|}_{i,j,k+1/2}^{n+1/2}\end{array}\right.$

$\left\{\begin{array}{c}{E}_x{|}_{i+1/2,j,k}^{n+1}=\frac{2-\mathrm{\Δt}\·b_z{|}_k}{2+\mathrm{\Δt}\·b_z{|}_k}{E}_x{|}_{i+1/2,j,k}^n+\frac{2+\mathrm{\Δt}\·b_x{|}_{i+1/2}}{2+\mathrm{\Δt}\·b_z{|}_k}{P}_x{|}_{i+1/2,j,k}^{n+1}-\frac{2-\mathrm{\Δt}\·b_x{|}_{i+1/2}}{2+\mathrm{\Δt}\·b_z{|}_k}{P}_x{|}_{i+1/2,j,k}^n\\ E_y{|}_{i,j+1/2,k}^{n+1}=\frac{2-\mathrm{\Δt}\·b_x{|}_i}{2+\mathrm{\Δt}\·b_x{|}_i}{E}_y{|}_{i,j+1/2,k}^n+\frac{2+\mathrm{\Δt}\·b_y{|}_{j+1/2}}{2+\mathrm{\Δt}\·b_x{|}_i}{P}_y{|}_{i,j+1/2,k}^{n+1}-\frac{2-\mathrm{\Δt}\·b_y{|}_{j+1/2}}{2+\mathrm{\Δt}\·b_x{|}_i}{P}_y{|}_{i,j+1/2,k}^n\\ E_z{|}_{i,j,k+1/2}^{n+1}=\frac{2-\mathrm{\Δt}\·b_y{|}_j}{2+\mathrm{\Δt}\·b_y{|}_j}{E}_z{|}_{i,j,k+1/2}^n+\frac{2+\mathrm{\Δt}\·b_z{|}_{k+1/2}}{2+\mathrm{\Δt}\·b_y{|}_j}{P}_z{|}_{i,j,k+1/2}^{n+1}-\frac{2-\mathrm{\Δt}\·b_z{|}_{k+1/2}}{2+\mathrm{\Δt}\·b_y{|}_j}{P}_z{|}_{i,j,k+1/2}^n\end{array}\right.$

E ^IncBe the incident field, P _x, P _y, P _zBe intermediate variable, b _x, b _y, b _zFor the absorption parameter of absorption layer can be determined by following formula:

$b_x{|}_i=\left\{\begin{array}{cc}{b}_x^{\max }\&CenterDot;{\left(\frac{n_{\mathrm{absorb}}-i}{n_{\mathrm{absorb}}}\right)}^3& i\&le;n_{\mathrm{absorb}}\\ 0& n_{\mathrm{absorb}}<i<N_x-n_{\mathrm{absorb}}\\ b_x^{\max }\&CenterDot;{\left(\frac{i-(N_x-n_{\mathrm{absorb}})}{n_{\mathrm{absorb}}}\right)}^3& i\&GreaterEqual;N_x-n_{\mathrm{absorb}}\end{array}\right.$

$b_y{|}_j=\left\{\begin{array}{cc}{b}_y^{\max }\&CenterDot;{\left(\frac{n_{\mathrm{absorb}}-j}{n_{\mathrm{absorb}}}\right)}^3& j\&le;n_{\mathrm{absorb}}\\ 0& n_{\mathrm{absorb}}<j<N_y-n_{\mathrm{absorb}}\\ b_y^{\max }\&CenterDot;{\left(\frac{j-(N_y-n_{\mathrm{absorb}})}{n_{\mathrm{absorb}}}\right)}^3& j\&GreaterEqual;N_y-n_{\mathrm{absorb}}\end{array}\right.$

$b_z{|}_k=\left\{\begin{array}{cc}{b}_z^{\max }\&CenterDot;{\left(\frac{n_{\mathrm{absorb}}-k}{n_{\mathrm{absorb}}}\right)}^3& k\&le;n_{\mathrm{absorb}}\\ 0& n_{\mathrm{absorb}}<k<N_z-n_{\mathrm{absorb}}\\ b_z^{\max }\&CenterDot;{\left(\frac{k-(N_z-n_{\mathrm{absorb}})}{n_{\mathrm{absorb}}}\right)}^3& k\&GreaterEqual;N_z-n_{\mathrm{absorb}}\end{array}\right.$

Wherein

C is the light velocity in the vacuum, chooses n in an embodiment of the present invention _Absorb=8;

The space derivative of each magnetic-field component that relates in the above-mentioned time iteration formula calculates with the pseudo-spectral method of the Fourier of this intranodal, promptly only operates at the magnetic field data of this computing node, comprises the overlapping region and the Non-overlapping Domain in magnetic field.Smooth function w (l), wherein 1≤l≤n be multiply by with the magnetic field of overlapping region in elder generation before pseudo-spectral method is asked the space derivative in magnetic field with Fourier _Overlap, the preceding half section value of w (l) is 1, and the value dullness of second half section drops to 0 (the qualitative Fig. 4 of opinion) from 1 glossily, and the magnetic field energy after requiring so to handle satisfies following three conditions: 1. the data in the Non-overlapping Domain are constant; 2. the data on this node satisfy periodic boundary condition, i.e. the data of Fourier transform operation about two end points be 0; 3. it is continuous that the first order derivative of magnetic field data is striden interphase.Repeatedly numerical experimentation shows that result's precision is insensitive to the precise forms of smooth function, only must satisfy that function upper end is smooth to be connected to 1, and the lower end is smooth to drop to 0.Choose smooth function in an embodiment of the present invention for (, establishing n for the programming convenience _OverlapBe even number)

$w\left(l\right)=\left\{\begin{array}{cc}1& 1\&le;l\&le;n_{\mathrm{overlap}}/2\\ \cos [\mathrm{\&pi;}(l-n_{\mathrm{overlap}}/2)/n_{\mathrm{overlap}}]& n_{\mathrm{overlap}}/2<l\&le;n_{\mathrm{overlap}}\end{array}\right.;$

L magnetic field data of overlapping region is multiplied by functional value w (l), makes the magnetic field data of overlapping region drop to 0 glossily from the end to end that is connected with Non-overlapping Domain.For example Fig. 5 is for asking the synoptic diagram of space derivative on the direction ξ (ξ can be x, y or z) in certain computing node, and arrow is represented the position in magnetic field, and round dot is represented the position of electric field, and the space derivative in magnetic field, electric field location place can be calculated as follows:

Further choose n in an embodiment of the present invention _Overlap=10, the form of smooth function is reduced to

$w\left(l\right)=\left\{\begin{array}{cc}1& 1\&le;l\&le;5\\ \cos [\mathrm{\&pi;}(l-5)/10]& 6\&le;l\&le;10\end{array}\right.$

Smooth function is affacted on the overlapping region in magnetic field, the magnetic field after the effect is

$\tilde{H}\left(i\right)=\left\{\begin{array}{cc}H\left(i\right)\&CenterDot;w({\mathrm{<figure-callout\; id="1"\; label="i"\; filenames="20101011377661000032.png,HSA00000037802600021.png"\; state="\{\{state\}\}">i</figure-callout>}}_1+10-i)& i_1\&le;i\&le;i_1+9\\ H\left(i\right)& i_1+10\&le;i\&le;i_2-10\\ H\left(i\right)\&CenterDot;w(i-{\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="20101011377661000032.png,HSA00000037802600031.png"\; state="\{\{state\}\}">i</figure-callout>}}_2+10)& i_2-9\&le;i\&le;i_2\end{array}\right.,$

Utilize expression to try to achieve the space derivative of magnetic field at the electric field location place

$\frac{\&PartialD;H}{\&PartialD;\mathrm{\&xi;}}=F^{-1}\left[{\mathrm{jk}}_{\mathrm{\&xi;}}{e}^{{\mathrm{jk}}_{\mathrm{\&xi;}}\mathrm{\&Delta;\&xi;}/2}F\left(\tilde{H}\right)\right],$

Wherein j is that imaginary unit is j ²=-1;

(5) each computing node receives the electric field data of the overlapping region that each neighborhood calculation node sends over, and sends the electric field data that the overlapping region of self Non-overlapping Domain and this neighborhood calculation node partially overlaps to each neighborhood calculation node;

(6) each computing node utilizes the time iteration formula in magnetic field to calculate in the Non-overlapping Domain next step magnetic field;

The time iteration formula in magnetic field is divided into two zones:

Non-absorption layer,

$\left\{\begin{array}{c}{H}_x{|}_{i,j+1/2,k+1/2}^{n+3/2}=H_x{|}_{i,j+1/2,k+1/2}^{n+1/2}+\frac{\mathrm{\&Delta;t}}{{\mathrm{\&mu;}}^{\&prime;}{|}_{i,j+1/2,k+1/2}}[(\frac{\&PartialD;E_y}{\&PartialD;z}-\frac{\&PartialD;E_z}{\&PartialD;y}){|}_{i,j+1/2,k+1/2}^{n+1}-({\mathrm{\&mu;}}^{\&prime;}{|}_{i,j+1/2,k+1/2}-{\mathrm{\&mu;}}_0^{\&prime;})\&CenterDot;\left(\frac{\&PartialD;H_{x,\mathrm{inc}}}{\&PartialD;t}\right){|}_{i,j+1/2,k+1/2}^{n+1}]\\ H_y{|}_{i+1/2,j,k+1/2}^{n+3/2}=H_y{|}_{i+1/2,j,k+1/2}^{n+1/2}+\frac{\mathrm{\&Delta;t}}{{\mathrm{\&mu;}}^{\&prime;}{|}_{i+1/2,j,k+1/2}}[(\frac{\&PartialD;E_z}{\&PartialD;x}-\frac{\&PartialD;E_x}{\&PartialD;z}){|}_{i+1/2,j,k+1/2}^{n+1}-({\mathrm{\&mu;}}^{\&prime;}{|}_{i+1/2,j,k+1/2}-{\mathrm{\&mu;}}_0^{\&prime;})\&CenterDot;\left(\frac{\&PartialD;H_{y,\mathrm{inc}}}{\&PartialD;t}\right){|}_{i+1/2,j,k+1/2}^{n+1}]\\ H_z{|}_{i+1/2,j+1/2,k}^{n+3/2}=H_z{|}_{i+1/2,j+1/2,k}^{n+1/2}+\frac{\mathrm{\&Delta;t}}{{\mathrm{\&mu;}}^{\&prime;}{|}_{i+1/2,j+1/2,k}}[(\frac{\&PartialD;E_x}{\&PartialD;y}-\frac{\&PartialD;E_y}{\&PartialD;x}){|}_{i+1/2,j+1/2,k}^{n+1}-({\mathrm{\&mu;}}^{\&prime;}{|}_{i+1/2,j+1/2,k}-{\mathrm{\&mu;}}_0^{\&prime;})\&CenterDot;\left(\frac{\&PartialD;H_{z,\mathrm{inc}}}{\&PartialD;t}\right){|}_{i+1/2,j+1/2,k}^{n+1}]\end{array}\right.$

Absorption layer,

$\left\{\begin{array}{c}{B}_x{|}_{i,j+1/2,k+1/2}^{n+3/2}=\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_y{|}_{j+1/2}}{2+\mathrm{\&Delta;t}\&CenterDot;b_y{|}_{j+1/2}}{B}_x{|}_{i,j+1/2,k+1/2}^{n+1/2}+\frac{2\mathrm{\&Delta;t}}{(2+\mathrm{\&Delta;t}\&CenterDot;b_y{|}_{j+1/2}){\mathrm{\&mu;}}_0^{\&prime;}}(\frac{\&PartialD;E_y}{\&PartialD;z}-\frac{\&PartialD;E_z}{\&PartialD;y}){|}_{i,j+1/2,k+1/2}^{n+1}\\ B_y{|}_{i+1/2,j,k+1/2}^{n+3/2}=\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_z{|}_{k+1/2}}{2+\mathrm{\&Delta;t}\&CenterDot;b_z{|}_{k+1/2}}{B}_y{|}_{i+1/2,j,k+1/2}^{n+1/2}+\frac{2\mathrm{\&Delta;t}}{(2+\mathrm{\&Delta;t}\&CenterDot;b_z{|}_{k+1/2}){\mathrm{\&mu;}}_0^{\&prime;}}(\frac{\&PartialD;E_z}{\&PartialD;x}-\frac{\&PartialD;E_x}{\&PartialD;z}){|}_{i+1/2,j,k+1/2}^{n+1}\\ B_z{|}_{i+1/2,j+1/2,k}^{n+3/2}=\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_x{|}_{i+1/2}}{2+\mathrm{\&Delta;t}\&CenterDot;b_x{|}_{i+1/2}}{B}_z{|}_{i+1/2,j+1/2,k}^{n+1/2}+\frac{2\mathrm{\&Delta;t}}{(2+\mathrm{\&Delta;t}\&CenterDot;b_x{|}_{i+1/2}){\mathrm{\&mu;}}_0^{\&prime;}}(\frac{\&PartialD;E_x}{\&PartialD;y}-\frac{\&PartialD;E_y}{\&PartialD;x}){|}_{i+1/2,j+1/2,k}^{n+1}\end{array}\right.$

$\left\{\begin{array}{c}{H}_x{|}_{i,j+1/2,k+1/2}^{n+3/2}=\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_z{|}_{k+1/2}}{2+\mathrm{\&Delta;t}\&CenterDot;b_z{|}_{k+1/2}}{H}_x{|}_{i,j+1/2,k+1/2}^{n+1/2}+\frac{2+\mathrm{\&Delta;t}\&CenterDot;b_x{|}_i}{2+\mathrm{\&Delta;t}\&CenterDot;b_z{|}_{k+1/2}}{B}_x{|}_{i,j+1/2,k+1/2}^{n+3/2}-\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_x{|}_i}{2+\mathrm{\&Delta;t}\&CenterDot;b_z{|}_{k+1/2}}{B}_x{|}_{i,j+1/2,k+1/2}^{n+1/2}\\ H_y{|}_{i+1/2,i,k+1/2}^{n+3/2}=\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_x{|}_{i+1/2}}{2+\mathrm{\&Delta;t}\&CenterDot;b_x{|}_{i+1/2}}{H}_y{|}_{i+1/2,j,k+1/2}^{n+1/2}+\frac{2+\mathrm{\&Delta;t}\&CenterDot;b_y{|}_j}{2+\mathrm{\&Delta;t}\&CenterDot;b_x{|}_{i+1/2}}{B}_y{|}_{i+1/2,j,k+1/2}^{n+3/2}-\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_y{|}_j}{2+\mathrm{\&Delta;t}\&CenterDot;b_x{|}_{i+1/2}}{B}_y{|}_{i+1/2,j,k+1/2}^{n+1/2}\\ H_z{|}_{i+1/2,j+1/2,k}^{n+3/2}=\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_y{|}_{j+1/2}}{2+\mathrm{\&Delta;t}\&CenterDot;b_y{|}_{j+1/2}}{H}_z{|}_{i+1/2,j+1/2,k}^{n+1/2}+\frac{2+\mathrm{\&Delta;t}\&CenterDot;b_z{|}_k}{2+\mathrm{\&Delta;t}\&CenterDot;b_y{|}_{j+1/2}}{B}_z{|}_{i+1/2,j+1/2,k}^{n+3/2}-\frac{2-\mathrm{\&Delta;t}\&CenterDot;b_z{|}_k}{2+\mathrm{\&Delta;t}\&CenterDot;b_y{|}_{j+1/2}}{B}_z{|}_{i+1/2,j+1/2,k}^{n+1/2}\end{array}\right.$

H ^IncFor the incident field distributes, B _x, B _y, B _zBe intermediate variable, b _x, b _y, b _zIdentical with parameter in above-mentioned (4);

The space derivative of each electric field component that relates in the above-mentioned time iteration formula calculates with the pseudo-spectral method of the Fourier of this intranodal, promptly only operates at the electric field data of this computing node, comprises the overlapping region and the Non-overlapping Domain of electric field.Earlier the electric field of overlapping region being multiply by smooth function before pseudo-spectral method is asked the electric field space derivative with Fourier, is the consistance that guarantees electric field and magnetic field space derivative calculations, adopt here with above-mentioned (4) in identical smooth function w (l).L electric field data of overlapping region are multiplied by functional value w (l), make the electric field data of overlapping region drop to 0 glossily from the end to end that is connected with Non-overlapping Domain.Magnetic field position place electric field space derivative can calculate as follows among Fig. 5:

Smooth function is affacted on the overlapping region of electric field, the electric field after the effect is

$\tilde{E}\left(i\right)=\left\{\begin{array}{cc}E\left(i\right)\&CenterDot;w({\mathrm{<figure-callout\; id="1"\; label="i"\; filenames="20101011377661000032.png,HSA00000037802600021.png"\; state="\{\{state\}\}">i</figure-callout>}}_1+10-i)& i_1\&le;i\&le;i_1+9\\ E\left(i\right)& i_1+10\&le;i\&le;i_2-10\\ E\left(i\right)\&CenterDot;w(i-{\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="20101011377661000032.png,HSA00000037802600031.png"\; state="\{\{state\}\}">i</figure-callout>}}_2+10)& i_2-9\&le;i\&le;i_2\end{array}\right.,$

Utilize expression to try to achieve the space derivative of electric field at the magnetic field position place

$\frac{\&PartialD;E}{\&PartialD;\mathrm{\&xi;}}=F^{-1}\left[{\mathrm{jk}}_{\mathrm{\&xi;}}{e}^{{\mathrm{jk}}_{\mathrm{\&xi;}}\mathrm{\&Delta;\&xi;}/2}F\left(\tilde{E}\right)\right];$

(7) receive the magnetic field data of the overlapping region that each neighborhood calculation node sends over, and send the magnetic field data that the overlapping region of self Non-overlapping Domain and this neighborhood calculation node partially overlaps to each neighborhood calculation node;

(8) judgement time goes on foot whether n is n ₁Or n ₂If:, then write down the electromagnetic field in the Non-overlapping Domain in each computing node this moment; If not, then enter step (9);

(9) judge whether n＜N sets up: if, n=n+1 and return step (4) then; If not, then enter step (10);

(10) all the other computing nodes except first computing node are with n ₁Step and n ₂The electromagnetic field that writes down during the step sends to first computing node; First computing node receives the data that send over from all the other each computing nodes, and distributes with the electromagnetic field that the data of himself record are merged in the whole virtual space, with n ₁The electromagnetic field in step space distributes and writes output file " result1 ", with n ₂The electromagnetic field in step space distributes and writes output file " result2 ";

(11) stop described MPI process.

Three, the scattered field that utilizes extrapolation far field, near field program to calculate the far field distributes;

1. read output file " result1 " and " result2 " of cluster computer, obtain in the space near field electromagnetic field at n ₁Step and n ₂Distribution during the step

With

The plural form of scattered field is obtained by following formula in the near field:

$\stackrel{\&RightArrow;}{E}={\stackrel{\&RightArrow;}{E}}_{n_1}-j{\stackrel{\&RightArrow;}{E}}_{{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="20101011377661000032.png,HSA00000037802600031.png"\; state="\{\{state\}\}">n</figure-callout>}}_2}$

$\stackrel{\&RightArrow;}{H}={\stackrel{\&RightArrow;}{H}}_{n_1}-j{\stackrel{\&RightArrow;}{H}}_{{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="20101011377661000032.png,HSA00000037802600031.png"\; state="\{\{state\}\}">n</figure-callout>}}_2};$

2. according to the requirement of extrapolation far field, the near field of standard program, a rectangular sealing surface is set (by six tangent plane x=x ₁, x=x ₂, y=y ₁, y=y ₂, z=z ₁, z=z ₂Surround, and can surround regional Ω);

3. interpolation obtains magnetic field vertical with this tangent plane normal direction on above-mentioned each tangent plane; As being acquisition face x=x ₁The magnetic-field component H of last y direction _y, at first to being positioned at face x=x ₁All H on the normal _yObtain its fourier spectrum as one dimensional fourier transform, then at the superior superior displacement factor of each spectrum component e ^{Jk Δ x/2}, k is corresponding discrete wave vector, then frequency spectrum is made the result of inverse-Fourier transform, promptly

$H_y^{\mathrm{sca}}{|}_{x=x_1}=F^{-1}\left(e^{\mathrm{jk\&Delta;x}/2}F\left(H_y^{\mathrm{sca}}\right)\right){|}_{x=x_1-\mathrm{\&Delta;x}/2};$

4. equivalent face electric current on the tangent plane and equivalent face magnetic current are defined as respectively

With

Wherein

For being unit normal vector outside on this tangent plane; Respectively each equivalent face current component on each tangent plane and each equivalent face magnetic current component are made two-dimension fourier transform, according to the frequency spectrum that obtains with they on the tangent plane of place the distribution expansion in series and form; For example at z=z ₁The component of equivalent surface current x direction is on the face

$J_x(x,y,z_1)=\underset{m,n}{\mathrm{\&Sigma;}}\widetilde{J_x}(m,n)e^{j(k_{x,m}x+k_{y,n}y)}$

For to whole z=z ₁Equivalent face electric current J on the tangent plane _xSpectrum component when making two-dimension fourier transform, k _{X, m}And k _{Y, n}Be respectively corresponding discrete wave vector;

5. for obtaining direction

The far field need to calculate earlier intermediate quantity

With Wherein

$\left\{\begin{array}{c}\stackrel{\&RightArrow;}{N}\left(\hat{r}\right)=\underset{s}{\&Integral;\&Integral;}(J_x\hat{x}+J_y\hat{y}+J_z\hat{z})e^{{\mathrm{jkr}}^{\&prime;}\cos \mathrm{\&psi;}}{\mathrm{ds}}^{\&prime;}\\ \stackrel{\&RightArrow;}{L}\left(\hat{r}\right)=\underset{s}{\&Integral;\&Integral;}(M_x\hat{x}+M_y\hat{y}+M_z\hat{z})e^{{\mathrm{jkr}}^{\&prime;}\cos \mathrm{\&psi;}}{\mathrm{ds}}^{\&prime;}\end{array}\right.,$

e ^{Jkr ' cos ψ}For to direction The relative phase delay amount at integration bin place during propagation, owing to equivalent face electric current and equivalent face magnetic current are launched into the form of clear and definite sum of series on each face, following formula can be resolved integration; For example to z=z ₁On the face

X durection component J _xIntegration is

${\&Integral;}_{x_1}^{x_2}{\&Integral;}_{y_1}^{y_2}{J}_x(x,y,z_1)e^{j(k_{x}x+k_{y}y+k_z{z}_1)}\mathrm{dxdy}$

$={\&Integral;}_{x_1}^{x_2}{\&Integral;}_{y_1}^{{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="20101011377661000032.png,HSA00000037802600031.png"\; state="\{\{state\}\}">y</figure-callout>}}_2}\underset{m,n}{\mathrm{\&Sigma;}}\widetilde{J_x}(m,n)e^{j(k_{x,m}x+k_{y,n}y)}{e}^{j(k_{x}x+k_{y}y+k_z{z}_1)}\mathrm{dxdy}$

$=e^{{\mathrm{<figure-callout\; id="1"\; label="jk\; z\; z"\; filenames="20101011377661000032.png,HSA00000037802600021.png"\; state="\{\{state\}\}">jk</figure-callout>}}_z{z}_{}{}_1}\underset{m,n}{\mathrm{\&Sigma;}}\left(\widetilde{J_x}(m,n)\frac{e^{j(k_{x,m}+k_x)x_2}-e^{j(k_{x,m}+k_x)x_1}}{j(k_{x,m}+k_x)}\frac{e^{j(k_{y,n}+k_y)y_2}-e^{j(k_{y,n}+k_y)y_1}}{j(k_{y,n}+k_y)}\right)$

K wherein _x, k _yAnd k _zBe respectively electromagnetic wave to direction

Wave vector during propagation is at the component of x, y and z direction.

6. calculate by following formula at last

The far field scattered field of direction

$\left\{\begin{array}{c}{E}_r=0\\ E_{\mathrm{\&theta;}}=-\frac{\mathrm{jk}{e}^{-\mathrm{<figure-callout\; id="4"\; label="jkr"\; filenames="20101011377661000032.png,HSA00000037802600041.png"\; state="\{\{state\}\}">jkr</figure-callout>}}}{4\mathrm{\&pi;r}}(L_{\mathrm{\&phi;}}+{\mathrm{\&eta;}}_0{N}_{\mathrm{\&theta;}})\\ E_{\mathrm{\&phi;}}=+\frac{{\mathrm{jke}}^{-\mathrm{<figure-callout\; id="4"\; label="jkr"\; filenames="20101011377661000032.png,HSA00000037802600041.png"\; state="\{\{state\}\}">jkr</figure-callout>}}}{4\mathrm{\&pi;r}}(L_{\mathrm{\&theta;}}-{\mathrm{\&eta;}}_0{N}_{\mathrm{\&phi;}})\end{array}\right.$

$\left\{\begin{array}{c}{H}_r=0\\ H_{\mathrm{\&theta;}}=+\frac{{\mathrm{jke}}^{-\mathrm{<figure-callout\; id="4"\; label="jkr"\; filenames="20101011377661000032.png,HSA00000037802600041.png"\; state="\{\{state\}\}">jkr</figure-callout>}}}{4\mathrm{\&pi;r}}(N_{\mathrm{\&phi;}}-\frac{L_{\mathrm{\&theta;}}}{{\mathrm{\&eta;}}_0})\\ H_{\mathrm{\&phi;}}=-\frac{{\mathrm{jke}}^{-\mathrm{<figure-callout\; id="4"\; label="jkr"\; filenames="20101011377661000032.png,HSA00000037802600041.png"\; state="\{\{state\}\}">jkr</figure-callout>}}}{4\mathrm{\&pi;r}}(N_{\mathrm{\&theta;}}+\frac{L_{\mathrm{\&phi;}}}{{\mathrm{\&eta;}}_0})\end{array}\right..$

Be the thought that the pseudo-spectral method of parallel processing electromagnetic field time domain is decomposed in further explaination overlapping region of the present invention, in Fig. 6, provided the synoptic diagram of overlapping region decomposition method, explain with the example that is decomposed in ξ direction zone.Shown in Fig. 6-1, the interval of zone on the ξ axle is [a, k], and width is N _ξNow will be divided into 4 computing nodes (computing node A, computing node B, computing node C and computing node D) of cluster computer, the width of data is L in each node _ξ, overlapping region thickness is n _OverlapShown in Fig. 6-2, at first [a, k] is divided into 10 sections [a, b], [b, c], [c, d], [d, e], [e, f], [f, g], [g, h], [h, i], [i, j] and [j, k], each section width is respectively L _ξ-2n _Overlap, n _Overlap, n _Overlap, L _ξ-4n _Overlap, n _Overlap, n _Overlap, L _ξ-4n _Overlap, n _Overlap, n _OverlapAnd L _ξ-2n _OverlapShown in Fig. 6-3, zone [a, d], [b, g], [e, j] and [h, k] are dispensed to computing node A, computing node B, computing node C and computing node D successively then.In computing node A, the electromagnetic field in zone [a, c] calculates with the pseudo-spectral method of electromagnetic field time domain, and zone [c, d] receives the corresponding segment data that is sent by computing node B; In computing node B, the electromagnetic field in zone [c, f] calculates with the pseudo-spectral method of electromagnetic field time domain, and zone [b, c] receives the corresponding segment data that is sent by computing node A, and zone [f, g] receives the corresponding segment data that is sent by computing node C; In computing node C, the electromagnetic field in zone [f, i] calculates with the pseudo-spectral method of electromagnetic field time domain, and zone [e, f] receives the corresponding segment data that computing node B sends, and zone [i, j] receives the corresponding segment data that computing node D sends; In computing node D, the electromagnetic field in zone [i, k] calculates with the pseudo-spectral method of electromagnetic field time domain, and zone [h, i] receives the corresponding segment data that computing node C sends.Fig. 6-4 has provided when asking the electromagnetic field space derivative in each computing node the synoptic diagram of the factor of taking advantage of on its data successively: factor perseverance is 1 in non-overlapped district, and electromagnetic field data is constant; Outwards reduce to 0 smoothly in the overlapping region factor by Non-overlapping Domain and border, overlapping region, the electromagnetic field data after the effect also drops to 0 smoothly.

Embodiment 1:

The polarization characteristic of rear orientation light during present embodiment research polarized light incident biological tissue model, thus parameter provided for the polarization gating technology in the optical detection biological tissue.Biological tissue's model adopts discrete scatterer model, as shown in Figure 7, and promptly at certain space volume range (l _x* l _y* l _z) the uniform medium microsphere of interior stochastic distribution.The radius of microballoon is 1 μ m, and refractive index is 1.59; Microballoon is a water with the external space, and refractive index is 1.33.Incident light is that its wavelength in free space is 785nm along the plane wave of the x directional ray polarization of Z direction input.

Common mode has been intended the organize models of four groups of different sizes in the present embodiment, and the number of microballoon in the bulk of each group model, the model, data length is listed in the table below in the width of grid and each computing node when cluster computer for simulating:

Model	The microballoon number	??l x??(μm)	??l y??(μm)	??l z??(μm)	??N x	??N y	??N z	??L x	??L y	??L z
											??1	??1735	??100	??100	??25	??1092	??1092	??300	??288	??288	??160
??2	??3470	??100	??100	??50	??1092	??1092	??556	??288	??288	??288
											??3	??5205	??100	??100	??75	??1092	??1092	??876	??288	??288	??448
??4	??6940	??100	??100	??100	??1092	??1092	??1132	??288	??288	??576

Chosen identical grid cell length (Δ x=Δ y=Δ z=0.098 μ m) in this external each group model, all simulations utilize 32 computing nodes of cluster computer, form 4 * 4 * 2 three-dimensional topology structure.

Grid comprises biological tissue's model of innermost layer, outermost absorption layer and between the two middle layer in the present embodiment 1.The thickness scalable in middle layer is so that parameter L _x, L _yAnd L _zSize meet the characteristic length of fast fourier transformation algorithm.

Fig. 8 is the process flow diagram of present embodiment 1:

(1) i=1 is set;

(2) j=1 is set;

(3) generate j realization of i group model: in the spatial dimension of i group model qualification, produce the microballoon position at random, guarantee all microballoon non-overlapping copies and in the space, be evenly distributed;

(4) set up j of i group model model and the generation model file of realizing;

(5) model that obtains in the step (4) is inputed to each computing node of cluster computer; The propagation condition of the pseudo-spectral method simulation plane wave of the electromagnetic field time domain that the cluster computer utilization has improved in this biological tissue's model; The data file that electromagnetic field in the cluster computer output record space behind the electromagnetic field distributional stability distributes;

(6) data file that obtains of read step (5) calculates the far field scattered field with extrapolation far field, near field program and distributes

With

(7) calculate scattering angle respectively in 175 °≤θ≤180 °, the total light intensity and with incident light polarization direction vertical total light intensity parallel with the incident light polarization direction in 0 °≤φ≤360 ° subtracted each other both and obtained polarization differential signal Δ I _{I, j}:

The electric field of the back scattering field that obtains in the step (6) is resolved into and parallel with the vertical two parts of the electric field of incident field, and resolve into and parallel with the vertical two parts in the magnetic field of incident field in magnetic field

$\left\{\begin{array}{c}{\stackrel{\&RightArrow;}{E}}_{i,j}(\mathrm{\&theta;},\mathrm{\&phi;})={\stackrel{\&RightArrow;}{E}}_{\left|\right|,i,j}(\mathrm{\&theta;},\mathrm{\&phi;})+{\stackrel{\&RightArrow;}{E}}_{\&perp;,i,j}(\mathrm{\&theta;},\mathrm{\&phi;})\\ {\stackrel{\&RightArrow;}{H}}_{i,j}(\mathrm{\&theta;},\mathrm{\&phi;})={\stackrel{\&RightArrow;}{H}}_{\left|\right|,i,j}(\mathrm{\&theta;},\mathrm{\&phi;})+{\stackrel{\&RightArrow;}{H}}_{\&perp;,i,j}(\mathrm{\&theta;},\mathrm{\&phi;})\end{array}\right.$

Thereby obtain

$\mathrm{\&Delta;}{I}_{i,j}=\frac{1}{2}\&Integral;\&Integral;\left(\right|{\stackrel{\&RightArrow;}{E}}_{\left|\right|,i,j}(\mathrm{\&theta;},\mathrm{\&phi;})\&times;{{\stackrel{\&RightArrow;}{H}}^{*}}_{\left|\right|,i,j}(\mathrm{\&theta;},\mathrm{\&phi;})|-|{\stackrel{\&RightArrow;}{E}}_{\&perp;,i,j}(\mathrm{\&theta;},\mathrm{\&phi;})\&times;{{\stackrel{\&RightArrow;}{H}}^{*}}_{\&perp;,i,j}(\mathrm{\&theta;},\mathrm{\&phi;})\left|\right)\sin \mathrm{\&theta;d\&theta;d\&phi;}.$

(8) judge whether j＜5 set up: if j=j+1 also returns step (3); If not, then enter step (9);

(9) the mean difference sub-signal of calculating i group model.

$\mathrm{\&Delta;}{I}_i=\frac{1}{5}\underset{j=1}{\overset{5}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{I}_{i,j}$

(10) judge whether i＜4 set up: if i=i+1 also returns step (2); If not, finish whole simulation.

Fig. 9 is the analog result of present embodiment, and horizontal ordinate is calculated by the Mie scattering theory for the optical thickness of each group biological tissue model correspondence, and ordinate is pairing model back scattering strength of differential signal.As can be seen from Figure 9, differential signal is saturated greater than 4 o'clock at optical thickness.Present embodiment has illustrated by detecting in the textura epidermoidea of several scattering lengths that the back scattering differential signal can be limited to sensing range biological tissue.

Embodiment 2:

Backscattering during present embodiment 2 analog computation light wave incident biological tissue models strengthens phenomenon.Biological tissue's model construction method in the present embodiment is consistent with biological tissue's model construction method among the embodiment 1, and each parameter is identical with the 2nd group model among the embodiment 1.Present embodiment is for suppressing the speckle of coherent light and biological tissue's model generation, and the result of the light incident of different frequency is done on average.

Figure 10 is the process flow diagram of present embodiment 2:

(1) generates organize models: in the spatial dimension that limits, produce the microballoon position at random, guarantee all microballoon non-overlapping copies and in the space, be evenly distributed;

(2) i=1 is set;

(3) calculate i the incident light frequency of simulating: f _i=(0.95+0.005 (i-1)) f ₀, f wherein ₀For the free space medium wavelength is the wavelength of 785nm light wave;

(4) setting up biological tissue's model is f in frequency _iDirectional light incident the time model and generation model file;

(5) model that obtains in the step (4) is inputed to each computing node of cluster computer; The propagation condition of the pseudo-spectral method simulation plane wave of the electromagnetic field time domain that the cluster computer utilization has improved in this biological tissue's model; The data file that electromagnetic field in the cluster computer output record space behind the electromagnetic field distributional stability distributes;

(6) data file that obtains of read step (5) calculates the back with extrapolation far field, near field program and distributes to the far field scattered field

With

(7) calculate the backscatter intensity I of different scattering angle according to the result of step (6) _i(θ)

$I_i\left(\mathrm{\&theta;}\right)=\frac{1}{2}\&Integral;|\stackrel{\&RightArrow;}{E_i}(\mathrm{\&theta;},\mathrm{\&phi;})\&times;{\stackrel{\&RightArrow;}{H_i}}^{*}(\mathrm{\&theta;},\mathrm{\&phi;})|\mathrm{d\&phi;};$

(8) judge whether i＜21 set up: if, i=i+1, and return step (3); If not, enter step (9);

The back scattering field is average when (7) getting different frequency incident

Figure 11 is the analog result of present embodiment 2: it is obvious that back scattering strengthens phenomenon, and obtained deriving from and organize the shallow-layer back scattering to strengthen angular range (～2 °).This simulation is applied in the optical detection biological tissue great importance is arranged back scattering being strengthened phenomenon.

It is worthy of note that the method for cluster computer for simulating electromagnetic wave propagation is not only applicable to light among above-mentioned two embodiment propagation simulation in biological tissue's model among the present invention, also be applicable to the propagation of any electromagnetic wave in complex space medium arbitrarily.

Claims (2)

1. the method for a cluster computer for simulating electromagnetic wave propagation is characterized in that, this method comprises the steps:

(1) the electromagnetic wave propagation problem of simulation is set up model and model of creation data file as required:

1. needing the electromagnetic wave propagation problem simulated for the electromagnetic wave of certain form incides behind the area of space Ω in the propagation in this zone, is certain single medium outside the regional Ω; Choose three mutually orthogonal directions arbitrarily

With

Satisfy

With Have the rectangular parallelepiped A of the minimum of an inclusion region Ω, each limit of A is parallel to

Or

Center with A is an initial point, With

The positive dirction that is respectively x axle, y axle and z axle is set up three-dimensional cartesian coordinate system;

2. choose grid cell length Δ x, Δ y and Δ z;

3. be 3D grid according to selected grid cell length with spatial division, be respectively N ' at the number of x, y and z direction grid cell _x, N ' _yAnd N ' _z, guarantee that rectangular parallelepiped A takes up space in 3D grid; Determine to meet following formula and meet and to make the fast fourier transformation algorithm minimum positive integer L of characteristic length efficiently _x, L _yAnd L _z:

L _x≥((N′ _x+2n _absorb+20)+2n _overlap(n _x-1))/n _x

L _y≥((N′ _y+2n _absorb+20)+2n _overlap(n _y-1))/n _y

L _z≥((N′ _z+2n _absorb+20)+2n _overlap(n _z-1))/n _z

N wherein _OverlapBe the thickness of overlapping region in the model, n _AbsorbBe the thickness of electromagnetic wave absorbing layer, n _x, n _yAnd n _zBe respectively the number of cluster computer three-dimensional topology structure at x, y and z direction calculating node; Calculate

N _x＝2(L _x-n _overlap)+(n _x-2)(l _x-2n _overlap)

N _y＝2(L _y-n _overlap)+(n _y-2)(L _y-2n _overlap)

n _z＝2(L _z-n _overlap)+(n _z-2)(L _z-2n _overap)；

4. with existing 3D grid at x axle positive dirction and the x axle negative direction n that stretches out respectively _X1And n _X2Layer grid cell obtains a new 3D grid, and the number of x, y and z direction grid cell is respectively N _x, N ' _yAnd N ' _z, wherein

n _X2=N _x-N ' _x-n _X1, function ceil (s) gets the minimum positive integer more than or equal to independent variable s;

5. with the 4. the new 3D grid that obtains of step at y axle positive dirction and the y axle negative direction n that stretches out respectively _Y1And n _Y2Layer grid cell obtains a new 3D grid, and the number of x, y and z direction grid cell is respectively N _x, N _yAnd N ' _z, wherein

n _Y2=N _y-N ' _y-n _Y1

6. with the 5. the new 3D grid that obtains of step at z axle positive dirction and the z axle negative direction n that stretches out respectively _Z1And n _Z2Layer grid cell obtains final 3D grid, and the number of x, y and z direction grid cell is respectively N _x, N _yAnd N _z, wherein n _Z2=N _z-N ' _z-n _Z1The occupied space of this 3D grid is the whole space of electromagnetic wave propagation simulation, outermost n _AbsorbLayer grid cell is electromagnetic wave absorbing layer;

7. the method according to the Yee grid disposes electric field and the position of magnetic field in each grid cell of 3D grid; Removal is positioned at all electromagnetic field components on positive dirction one side end face of 3D grid all directions, respectively corresponding three-dimensional matrice of each electromagnetic field component in the grid, and these matrixes are of similar shape: the length of x direction is N _x, the length of y direction is N _y, the length of z direction is N _z

8. the whole spatial division of utilizing the overlapping region decomposition method to simulate is piece, the number of piece equals the number of this direction calculating node in the cluster computer three-dimensional topology structure respectively on x, y and the z direction, the scope of each electromagnetic field component is determined by following method in each piece: the matrix that each electromagnetic field component is formed is divided into continuous submatrix respectively, the number of all directions submatrix equals the number of piece on this direction respectively, and the number of each submatrix element on direction ξ is L _ξ-n _OverlapOr L _ξ-2n _Overlap, when this submatrix is the head of ξ direction or odd amount in addition to the round number matrix, be L _ξ-n _Overlap, all the other situations are L _ξ-2n _Overlap, wherein ξ is x, y or z; The zone of each submatrix of gained is extended, principle is if this submatrix has in adjacent submatrix and this submatrix corresponding electromagnetic field component direction and this border tangent on certain border, then with this submatrix at this boundary coordinate of border vertical direction n that stretches out therewith _OverlapEach submatrix of above-mentioned each electromagnetic field component is one by one corresponding to each piece, and the zone before extending is the Non-overlapping Domain in this piece, and the zone of Yan Shening was the overlapping region afterwards;

9. access time, step delta t was

Wherein

T is the electromagnetic cycle, v _MaxThe maximum rate of in each medium of space, propagating for electromagnetic wave;

10. determine the time step n that iterations N and record electromagnetic field distribute ₁And n ₂

The specific inductive capacity and the magnetic permeability of each medium in the space are revised as ε ' respectively _h=α ε _hAnd μ ' _h=α μ _h, wherein

Be electromagnetic angular frequency, ε _hAnd μ _hBe respectively the specific inductive capacity and the magnetic permeability of h class medium reality, ε ' _hAnd μ ' _hBe respectively the h class medium specific inductive capacity and the magnetic permeability that in simulation, use;

The model of creation data file: described overlapping region decomposition method is cut apart each piece that obtains, set up the model data file respectively, and called after " model_i _b_ j _b_ k _b", i wherein _b, j _bAnd k _bBe respectively the sequence number of this piece, following content write the pairing model data file of each piece in x, y and z direction: in this piece each electromagnetic field component in the coordinate range of all directions overlapping region and Non-overlapping Domain, this piece in the Non-overlapping Domain each amended specific inductive capacity of electromagnetic field component position medium and magnetic permeability, grid cell at the length of all directions, time step, iterations, the time step n in the simulation ₁And n ₂, and the parameter of incident field;

(2) cluster computer carries out data processing:

Read the model data file, with the pseudo-spectral method Parallel Simulation of parallel electromagnetic field time domain electromagnetic wave propagation, output record time step n ₁And n ₂The destination file that the time space electromagnetic field distributes;

1. the model data file of resulting correspondence in each computing node difference read step (1) of cluster computer obtains simulating required parameter, and coordinate is (i in the cluster computer three-dimensional topology structure _p, j _p, k _p) the model data file of computing node correspondence be " model_i _p_ j _p_ k _p"; The space of each computing node each electromagnetic field component of memory allocated in internal memory;

2. utilize the electromagnetic wave propagation in the pseudo-spectral method virtual space of parallel electromagnetic field time domain;

3. export destination file " result1 " and " result2 ", this destination file " result1 " and " result2 " have write down n respectively ₁Step and n ₂The numerical value of step each electromagnetic field component of time space;

(3) scattered field that utilizes extrapolation far field, near field program to calculate the far field distributes:

1. destination file " result1 " and " result2 " of output in the read step (2) obtain the complex field of each electromagnetic field component in the near field;

2. in 3D grid, the sealing surface of a virtual enclosing region Ω is set, by six tangent plane x=x ₁, x=x ₂, y=y ₁, y=y ₂, z=z ₁And z=z ₂Surround;

3. interpolation obtains each vertical with this tangent plane normal direction on above-mentioned each tangent plane magnetic-field component;

4. on each tangent plane, respectively each equivalent face electric current and each equivalent face magnetic current are made two-dimension fourier transform, each the equivalent face electric current on each tangent plane and each equivalent face magnetic current are expanded into the form of sum of series;

5. resolve integration and obtain the direction of extrapolating

The far field time required intermediate quantity

With

6. by

With

Calculate direction

The far field scattered field With

2. the method for cluster computer for simulating electromagnetic wave propagation according to claim 1, it is characterized in that, the pseudo-spectral method of described parallel electromagnetic field time domain is by time iteration simulation electromagnetic wave propagation, and electromagnetic field is finished following four steps successively by step current time iteration each computing node to the process of next time step in the space:

(1) utilize the time iteration formula of electric field to ask the electric field of Non-overlapping Domain in next this computing node of time step; Wherein the space derivative of each magnetic-field component that relates in the time iteration formula of electric field calculates with the pseudo-spectral method of the Fourier of this intranodal, promptly only operate at the magnetic field data of this computing node, comprise its overlapping region and Non-overlapping Domain, and earlier on the magnetic field of overlapping region, be multiplied by a smooth function and make it drop to 0 smoothly from the end to end that is connected with Non-overlapping Domain;

(2) receive the electric field data of the overlapping region that each neighborhood calculation node sends over, and send the electric field data that the overlapping region of self Non-overlapping Domain and this neighborhood calculation node partially overlaps to each neighborhood calculation node;

(3) utilize the time iteration formula in magnetic field to ask the magnetic field of Non-overlapping Domain in next this computing node of time step; Wherein the space derivative of each electric field component that relates in the time iteration formula in magnetic field calculates with the pseudo-spectral method of the Fourier of this intranodal, promptly only operate at the electric field data of this computing node, comprise its overlapping region and Non-overlapping Domain, and earlier on the electric field of overlapping region, be multiplied by a smooth function and make it drop to 0 smoothly from the end to end that is connected with Non-overlapping Domain;

(4) receive the magnetic field data of the overlapping region that each neighborhood calculation node sends over, and send the magnetic field data that the overlapping region of self Non-overlapping Domain and this neighborhood calculation node partially overlaps to each neighborhood calculation node.

Images