JP2005024354A - Parallel sound field calculator - Google Patents

Parallel sound field calculator Download PDF

Info

Publication number
JP2005024354A
JP2005024354A JP2003189156A JP2003189156A JP2005024354A JP 2005024354 A JP2005024354 A JP 2005024354A JP 2003189156 A JP2003189156 A JP 2003189156A JP 2003189156 A JP2003189156 A JP 2003189156A JP 2005024354 A JP2005024354 A JP 2005024354A
Authority
JP
Japan
Prior art keywords
sound field
processing unit
calculation
sound
outputs
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
JP2003189156A
Other languages
Japanese (ja)
Inventor
Nobuo Makimoto
伸生 牧本
Kazuhiko Hino
一彦 日野
Kenichi Asano
健一 浅野
Hiroshi Takahashi
廣 高橋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hitachi Ltd
Original Assignee
Hitachi Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hitachi Ltd filed Critical Hitachi Ltd
Priority to JP2003189156A priority Critical patent/JP2005024354A/en
Publication of JP2005024354A publication Critical patent/JP2005024354A/en
Pending legal-status Critical Current

Links

Images

Abstract

<P>PROBLEM TO BE SOLVED: To simulate an underwater sound field at high speed as to the areas of remote sensing and underwater acoustic information processing. <P>SOLUTION: Correspondently to a sound wave propagation model, processing is divided for parallel execution. That is, a multitude of characteristic values and characteristic vectors are parallel calculated as to a normal mode method. As to an FDTD (Finite Difference Time Domain) method, a domain is divided and parallel calculation is performed while suitably exchanging boundary data. Sound field calculation can be speeded up to facilitate the actualization of repetitive simulation and real-time simulation. <P>COPYRIGHT: (C)2005,JPO&NCIPI

Description

【０００１】
【発明の属する技術分野】
本発明はリモートセンシングや水中音響情報処理の分野に関する。
【０００２】
【従来の技術】
各種の音響モデルを説明するため、まず水中音響で用いられる波動の方程式について説明する。座標系を図５のとおりとし、記号はｐを音圧、ρを密度、ｃを音速、ｔを時刻、ｋを波数とする。次の式は通常の波動方程式である。
【０００３】
【数１】

Figure 2005024354
水中音響の分野では、これに音源の定常性と等方性を仮定して極座標表示した次の方程式（Ｂｅｓｓｅｌ の方程式）がよく用いられる。
【０００４】
【数２】
Figure 2005024354
次に各種の音響モデルの概要を説明する。これらについては「海洋音響学会誌 Ｖｏｌ．２８ Ｎｏ．４」（平成１３年１０月）により詳しい説明がある。
【０００５】
ノーマルモード法とは、Ｂｅｓｓｅｌ の方程式において
【０００６】
【数３】
Figure 2005024354
のように深度と距離とを変数分離して
【０００７】
【数４】
Figure 2005024354
の形の固有値問題に帰着する方法である。数値的には、上記の固有値問題を差分化した行列の固有値・固有ベクトルを多数計算し、上記の変数分離式のように重ね合わせることで音場ｐ＝ｐ（ｒ，ｚ）を算出する。
【０００８】
ＰＥ（Ｐａｒａｂｏｌｉｃ Ｅｑｕａｔｉｏｎ； 放物型方程式）法とは、Ｂｅｓｓｅｌ の方程式を
【０００９】
【数５】
Figure 2005024354
のように Ｈａｎｋｅｌ 関数Ｈ０と摂動項Ψとに変数分離して変形し（ｋ０ は典型的な波数をあらわす定数）、
【００１０】
【数６】
Figure 2005024354
の方程式に帰着するものである。Ｘはｚのみについての２階の微分作用素である。ＰＥ法の名前は、上記方程式がｒを時間に見立てると熱方程式（放物型方程式）に近い形であることに由来する。Ｈａｎｋｅｌ 関数とは、Ｂｅｓｓｅｌ の方程式の解としてよく知られたものである。これについては日本数学会編「数学事典」岩波書店に説明がある。数値的には、上記方程式を差分化した３重対角行列係数の１次方程式を解き、ｒが小さいほうから大きいほうへ摂動Ψ＝Ψ（ｒ，ｚ）を順次算出する。
【００１１】
ＦＤＴＤ（Ｆｉｎｉｔｅ Ｄｉｆｆｅｒｅｎｃｅ Ｔｉｍｅ Ｄｏｍａｉｎ； 有限差分時間領域）法は、波動方程式を導くもととなる媒質の状態方程式
【００１２】
【数７】
Figure 2005024354
と媒質の運動方程式
【００１３】
【数８】
Figure 2005024354
とを連立させて解くものである。前者は音圧が変位に比例する Ｈｏｏｋ の法則をあらわす状態方程式、後者は媒質の加速度が音圧勾配に比例し粘性抵抗の分だけ減るすることを現す運動方程式である。ｐ＝ｐ（ｒ，ｚ，ｔ）は音圧、ν＝ν（ｒ，ｚ，ｔ）は媒質の速度、Ｋとηは別途定まる係数である。媒質速度νが決まると状態方程式により次の時刻の音圧ｐが決まり、音圧ｐ決まると運動方程式により次の時刻の媒質速度νが決まる構造になっている。数値的には、これらを差分化した ｌｅａｐ ｆｌｏｇ（蛙飛び）アルゴリズムと呼ばれるｐ，νの交互計算により音場ｐを時々刻々算出する。ノーマルモード法およびＰＥ法では時間非依存の音場ｐ＝ｐ（ｒ，ｚ）を算出するのに対し、ＦＤＴＤ法では時間依存の音場ｐ＝ｐ（ｒ，ｚ，ｔ）を算出する。
【００１４】
従来、電波伝播を対象にしたＦＤＴＤ法を並列化した例が打矢・柏「並列型スーパーコンピュータを用いたＦＤＴＤ並列計算」ＩＥＩＣＥ Ｔｒａｎｓ． ｖｏｌ．Ｊ８４．Ｃ Ｎｏ．１１ ２００１ に記載されている。しかし基本方程式の異なる音波伝播を対象にしたＦＤＴＤ法を並列化した例は見当たらない。またノーマルモード法やＰＥ法を並列化した例は見当たらない。
【００１５】
【非特許文献１】
「海洋音響学会誌 Ｖｏｌ．２８ Ｎｏ．４」（平成１３年１０月）
【非特許文献２】
日本数学会編「数学事典」岩波書店
【非特許文献３】
打矢・柏「並列型スーパーコンピュータを用いたＦＤＴＤ並列計算」ＩＥＩＣＥ Ｔｒａｎｓ． ｖｏｌ．Ｊ８４．Ｃ Ｎｏ．１１ ２００１
【００１６】
【発明が解決しようとする課題】
本発明は水中音響シミュレーションを高速におこなうことを課題とする。
【００１７】
【課題を解決するための手段】
ノーマルモード法については、多数の固有値・固有ベクトルを並列計算する。ＦＤＴＤ法については、領域を分割して適宜境界データを交換しながら並列計算する。
【００１８】
【発明の実施の形態】
第１の例として、ＦＤＴＤ法を並列化する装置を説明する。図１がその構成およびデータの流れをあらわしている。音源データ１１１は音源の周波数や位置などからなる。海洋１１２は音速プロファイルや吸収率プロファイルからなる。格子化の条件１１３は、解析の空間的および時間的範囲とその刻み幅、出力の空間的および時間的範囲とその刻み幅（通常は後者のほうがかなり粗い）などからなる。初期化処理部１０１では、海洋１１２と格子１１３の条件から定まる差分スキーム１１５と、格子から決まる初期音場ｐ・初期媒質速度場νの格子データ１１６とを出力する。ｐ、νの初期値は通常至る所０である。分割処理部１０２は、解析および出力の領域を分割し、また分割した解析領域を計算機に割り振る。ここで差分スキーム１１５は分割された領域に対応する分割された差分スキーム１１７に、同じく初期音場・初期媒質速度場１１６は分割された初期音場・媒質速度場１１８に、同じく解析用・出力用の格子１１６は分割された格子１１９に変換され、また分割された領域で音場計算をするためのＦＤＴＤサーバ１０６を起動する各計算機やそれらサーバ間の接続方法を記述するスクリプト１２０が生成・実行される。配布処理部では、分割された差分スキーム１１７、初期データ１１８、格子１１９を対応するＦＤＴＤサーバ１０６に配布する。部分音場計算処理部１０４内では、各ＦＤＴＤサーバ１０６が時々刻々の音場を計算し、隣接する部分領域から境界のデータを交換しながら全体領域の音場をもれなく計算する。回収統合処理部１０５では、上記各ＦＤＴＤサーバ１０６の音場計算結果を定期的に回収し、それらをつなげて１個の音場データに統合する。算出された音場１２１は必要に応じて加工・保存される。
【００１９】
境界データの交換の仕方を、図２をもちいて説明する。ある時刻ｔ（ｋ）において、各部分計算領域で正しい音場ｐ・媒質速度場νを保持しているとする。次の境界データ交換時刻ｔ（ｋ＋１）＝ｔ（ｋ）＋Δｔまでの処理にうつる前に、各部分計算領域は境界の内側の正しいデータを隣接領域（辺を接するもの）に提供し、同じく境界の外側の正しいデータを隣接領域から取得・保持する。音速をｃとすると、次のデータ交換時刻まで境界の内側で正しいデータを保持するには、境界からｃΔｔの厚さに相当する部分のデータを交換すればよい。ここで、まず上下方向のデータ交換をおこない、それが完了してから左右方向のデータ交換をおこなうようにする。こうすることで、図２の濃い色の部分のような、直接辺を接する領域が正しい値を保持していない部分についても、対角にある領域から正しい値を取得できる。なお境界データ交換の対象の領域の厚さ（符号２３０）〜ｃΔｔを大きくとるとデータ交換の回数が少なくて済むので、通信のオーバーヘッドが大きい場合に好適である。
【００２０】
第２の例として、ノーマルモード法を並列化する方法を図４をもちいて説明する。初期化処理４０１で差分スキームを決定したあとは、従来の技術の節でも説明したように、ノーマルモード法の計算の中心は３重対角行列の固有値・固有ベクトルを（小さいほうから）多数算出する部分４０２である。この部分を並列に実行し、結果を回収統合処理部４０３で集めて重ね合わせ、最終的に音場４１５を得る。したがってこの固有値・固有ベクトル４０２を互いに独立に算出できれば十分であり、以下ｋ番目の固有値・固有ベクトルを算出する方法の例を説明する。
【００２１】
対称行列Ａのλ’以下の固有値の数をＳ（Ａ，λ’）という記号であらわす。３重対角行列の場合、Ｓ（Ａ，λ’）の値は Ｓｔｕｒｍ の定理を用いて容易に計算できる。これについては 名取「数値解析とその応用」 １９９０ コロナ社 に説明があるとおりである。明らかにＳ（Ａ，−｜Ａ｜−ε）＝０、Ｓ（Ａ，｜Ａ｜＋ε）＝ｎ（次元）であり、λａ ＝ −｜Ａ｜−ε、λｂ ＝ ｜Ａ｜＋εはｋ番目の固有値の下界および上界になっている。そこでつぎの二分法の処理を繰り返す。
【００２２】
【数９】
Figure 2005024354
λａとλｂ とが他の固有値に比べｋ番目の固有値に十分近づいたら、つぎの逆反復法の処理を収束まで繰り返す。
【００２３】
【数１０】
Figure 2005024354
これにより、ｋ番目の固有値が縮退していない場合はｋ番目の固有値・固有ベクトルが得られる。なお固有値が縮退している場合は、すでに求まっている固有ベクトルで所定の初期値を Ｓｃｈｍｉｄｔ 直交化したものを逆反復法の初期値として用いれば、並列化の効果は損なわれるものの、独立な固有値・固有ベクトルを順次得ることができる。
【００２４】
第３の例として、第１の例の並列化ＦＤＴＤを中断・再開する方法を、図３、図１を用いて説明する。まず、図１の回収統合処理部１０５が中断要求を常時監視３１０するようにする。シグナルやスレッドを用いれば、ＦＤＴＤの計算を妨げずに中断要求を監視することができる。中断要求を検出した場合、部分計算同期処理３２０にうつる。ここではまず回収統合処理部１０５が各部分計算処理部１０６に中断要求を転送する。ｉ番目の部分計算処理部１０６は、この中断要求を検出した場合、ＦＤＴＤの計算をきりのよいところまで進めてから一時停止し、処理済のシミュレーション時刻ｔｉを回収統合処理部１０５に報告する。回収統合処理部１０５はこれらのシミュレーション時刻の最大値ｔｍａｘを計算し、各部分計算処理部１０６に通知する。各部分計算処理部１０６は、シミュレーション時刻ｔｍａｘまでＦＤＴＤの処理を通常どおりに進める。ここで音場回収統合処理３３０にうつる。各部分計算処理部１０６は、時刻ｔｍａｘにおいて保持している部分計算領域の音場（スカラー場）ｐおよび媒質速度場（ベクトル場）νを回収統合処理部１０５に出力する。回収統合処理部１０５は各部分計算処理部１０６からの上記データを回収・統合し、時刻ｔｍａｘにおける全体の音場ｐおよび媒質速度場νを作成する。ここで音場出力処理３４０にうつる。ここでは、上記の音場ｐ、媒質速度場ν、時刻ｔｍａｘの各データを所定のファイルに保存する。これで中断処理は完了である。再開の場合は、保存した音場ｐ、媒質速度場ν、時刻ｔｍａｘの各データを初期化処理部１０１で読み込めばよい。
【００２５】
【発明の効果】
以上のように、本発明により各種音響モデルにもとづき音場を高速に計算できる。
【図面の簡単な説明】
【図１】並列ＦＤＴＤ法による音場計算装置の構成図。
【図２】並列ＦＤＴＤ法による部分計算処理部間の境界データ交換の説明図。
【図３】並列ＦＤＴＤ法による音場計算装置の中断処理の説明図。
【図４】並列ノーマルモード法による音場計算装置の構成図。
【図５】音響モデルを記述する座標系を示す図。
【符号の説明】
１００…並列ＦＤＴＤ法による音場計算装置の構成およびデータの流れ、１０１…初期化処理部、１０２…分割処理部、１０３…配布処理部、１０４…部分計算処理部の集合、１０５…回収統合処理部、１０６…部分計算処理を担当する個々のＦＤＴＤサーバ、１１１…音源を規定するデータ、１１２…海洋などの媒質の場を規定するデータ（音速プロファイルや深度プロファイルなど）、１１３…差分化用の格子（範囲と刻み幅など）、１１４…トポロジー（部分計算における領域分割方法、部分計算領域と計算機とのマッピングなど）、１１５…差分スキーム（差分化方程式の係数）、１１６…初期音場（および初期速度場）、１１７…部分計算領域ごとに分割された差分スキーム、１１８…部分計算領域ごとに分割された初期音場（および初期速度場）、１１９…部分計算領域ごとに分割された差分化用格子、１２０…部分計算処理および回収統合処理の起動データ、１２１…全体の音場、２００…境界データ交換の説明、２１０…部分計算領域、２２０…データ交換領域、２３０…データ交換領域の厚さ、３００…音場計算中断の処理、３１０…中断要求検出処理、３２０…部分計算領域同期処理、３３０…音場回収統合処理、３４０…音場出力処理、３５０…中断時の音場データ、４００…並列ノーマルモード法による音場計算装置の構成およびデータの流れ、４０１…初期化処理部、４０２…固有値・固有ベクトル算出処理部、４０３…回収統合処理部、４１１…音源を規定するデータ、４１２…海洋などの媒質の場を規定するデータ（音速プロファイルなど）、４１３…差分化用の格子（範囲と刻み幅など）、４１４…差分スキーム（差分化方程式の係数）、４１５…音場。[0001]
BACKGROUND OF THE INVENTION
The present invention relates to the fields of remote sensing and underwater acoustic information processing.
[0002]
[Prior art]
In order to explain various acoustic models, first, equations of waves used in underwater acoustics will be described. The coordinate system is as shown in FIG. 5. Symbols are p for sound pressure, ρ for density, c for sound velocity, t for time, and k for wave number. The following equation is the usual wave equation.
[0003]
[Expression 1]
Figure 2005024354
In the field of underwater acoustics, the following equation (Bessel equation) expressed in polar coordinates on the assumption that the sound source is stationary and isotropic is often used.
[0004]
[Expression 2]
Figure 2005024354
Next, the outline of various acoustic models will be described. These are described in more detail in “Journal of Ocean Acoustics Vol.28 No.4” (October 2001).
[0005]
The normal mode method is Bessel's equation.
[Equation 3]
Figure 2005024354
Separate the variable of depth and distance as
[Expression 4]
Figure 2005024354
Is a method that reduces to an eigenvalue problem of the form Numerically, the sound field p = p (r, z) is calculated by calculating a large number of eigenvalues / eigenvectors of a matrix obtained by differentiating the above eigenvalue problem and superposing them as in the above variable separation formula.
[0008]
The PE (parabolic equation) method is the Bessel equation.
[Equation 5]
Figure 2005024354
As shown, the Hankel function H0 and the perturbation term Ψ are separated into variables (k0 is a constant representing a typical wavenumber)
[0010]
[Formula 6]
Figure 2005024354
This results in X is a second-order differential operator for z only. The name of the PE method is derived from the fact that the above equation is close to the heat equation (parabolic equation) when r is regarded as time. The Hankel function is well known as the solution of the Bessel equation. This is explained in the “Mathematical Encyclopedia” Iwanami Shoten edited by the Mathematical Society of Japan. Numerically, a linear equation of tridiagonal matrix coefficients obtained by differentiating the above equation is solved, and perturbation Ψ = Ψ (r, z) is sequentially calculated from the smaller r to the larger.
[0011]
The FDTD (Finite Difference Time Domain) method is a state equation of a medium from which a wave equation is derived.
[Expression 7]
Figure 2005024354
And the equation of motion of the medium
[Equation 8]
Figure 2005024354
And solve them together. The former is the equation of state that represents Hook's law in which the sound pressure is proportional to the displacement, and the latter is the equation of motion that indicates that the acceleration of the medium is proportional to the sound pressure gradient and decreases by the viscous resistance. p = p (r, z, t) is the sound pressure, ν = ν (r, z, t) is the velocity of the medium, and K and η are separately determined coefficients. When the medium speed ν is determined, the sound pressure p at the next time is determined by the state equation, and when the sound pressure p is determined, the medium speed ν at the next time is determined by the equation of motion. Numerically, the sound field p is calculated from time to time by alternating calculation of p and ν, which is called a “leap flag” algorithm in which these are differentiated. In the normal mode method and the PE method, a time-independent sound field p = p (r, z) is calculated, whereas in the FDTD method, a time-dependent sound field p = p (r, z, t) is calculated.
[0014]
Conventionally, an example of parallelization of the FDTD method for radio wave propagation has been described by Uchiya and Tsuji, “FDTD Parallel Computation Using a Parallel Supercomputer,” IEICE Trans. vol. J84. C No. 11 2001. However, there is no example of parallelizing the FDTD method for sound wave propagation with different basic equations. There are no examples of parallelizing the normal mode method or PE method.
[0015]
[Non-Patent Document 1]
"Journal of Ocean Acoustics Society Vol.28 No.4" (October 2001)
[Non-Patent Document 2]
Mathematical encyclopedia edited by the Japanese Mathematical Society, Iwanami Shoten [Non-patent Document 3]
Uchiya, Satoshi "FDTD Parallel Computation Using Parallel Supercomputer" IEICE Trans. vol. J84. C No. 11 2001
[0016]
[Problems to be solved by the invention]
An object of the present invention is to perform underwater acoustic simulation at high speed.
[0017]
[Means for Solving the Problems]
For the normal mode method, a large number of eigenvalues and eigenvectors are calculated in parallel. For the FDTD method, parallel calculation is performed while dividing the region and exchanging boundary data as appropriate.
[0018]
DETAILED DESCRIPTION OF THE INVENTION
As a first example, an apparatus for parallelizing the FDTD method will be described. FIG. 1 shows the configuration and the flow of data. The sound source data 111 includes the frequency and position of the sound source. The ocean 112 includes a sound velocity profile and an absorptivity profile. The grid condition 113 includes a spatial and temporal range of analysis and its step size, a spatial and temporal range of output, and its step size (usually the latter is considerably coarser). The initialization processing unit 101 outputs a difference scheme 115 determined from the conditions of the ocean 112 and the lattice 113 and lattice data 116 of the initial sound field p and the initial medium velocity field ν determined from the lattice. The initial values of p and ν are usually zero everywhere. The division processing unit 102 divides analysis and output areas, and allocates the divided analysis areas to computers. Here, the difference scheme 115 is divided into the divided difference scheme 117 corresponding to the divided area, and the initial sound field / initial medium velocity field 116 is also divided into the initial sound field / medium velocity field 118, and is also used for analysis / output. The grid 116 is converted into a divided grid 119, and a script 120 describing each computer that starts the FDTD server 106 for calculating the sound field in the divided area and a connection method between the servers is generated and generated. Executed. The distribution processing unit distributes the divided difference scheme 117, initial data 118, and grid 119 to the corresponding FDTD server 106. In the partial sound field calculation processing unit 104, each FDTD server 106 calculates a sound field from moment to moment, and calculates the sound field of the entire region while exchanging boundary data from adjacent partial regions. The collection integration processing unit 105 periodically collects the sound field calculation results of each of the FDTD servers 106, connects them, and integrates them into one sound field data. The calculated sound field 121 is processed and stored as necessary.
[0019]
A method of exchanging the boundary data will be described with reference to FIG. It is assumed that a correct sound field p and medium velocity field ν are held in each partial calculation region at a certain time t (k). Before going to the next boundary data exchange time t (k + 1) = t (k) + Δt, each partial calculation area provides the correct data inside the boundary to the adjacent area (the one that touches the edge), and also the boundary The correct data outside is acquired and retained from the adjacent area. Assuming that the speed of sound is c, in order to hold correct data inside the boundary until the next data exchange time, the data corresponding to the thickness of cΔt from the boundary may be exchanged. Here, data exchange in the vertical direction is first performed, and data exchange in the horizontal direction is performed after the data exchange is completed. By doing so, the correct value can be obtained from the diagonal region even for the portion where the region directly touching the side, such as the dark colored portion in FIG. 2, does not hold the correct value. Note that if the thickness (reference numeral 230) to cΔt of the region to be subjected to boundary data exchange is increased, the number of data exchanges can be reduced, which is preferable when communication overhead is large.
[0020]
As a second example, a method of parallelizing the normal mode method will be described with reference to FIG. After the difference scheme is determined in the initialization process 401, as described in the section of the prior art, the center of the calculation of the normal mode method is to calculate a large number of eigenvalues / eigenvectors of a tridiagonal matrix (from the smallest). Part 402. This part is executed in parallel, and the results are collected and overlapped by the collection and integration processing unit 403 to finally obtain a sound field 415. Therefore, it is sufficient if the eigenvalue / eigenvector 402 can be calculated independently of each other. An example of a method for calculating the kth eigenvalue / eigenvector will be described below.
[0021]
The number of eigenvalues less than or equal to λ ′ of the symmetric matrix A is represented by the symbol S (A, λ ′). In the case of a tridiagonal matrix, the value of S (A, λ ′) can be easily calculated using the theorem of Sturm. This is as explained in Natori "Numerical analysis and its application" 1990 Corona. Clearly, S (A, − | A | −ε) = 0, S (A, | A | + ε) = n (dimensions), and λa = − | A | −ε and λb = | A | + ε are k The lower and upper bounds of the eigenvalue. Therefore, the following bisection method is repeated.
[0022]
[Equation 9]
Figure 2005024354
When λa and λb are sufficiently close to the k-th eigenvalue compared to the other eigenvalues, the next inverse iteration process is repeated until convergence.
[0023]
[Expression 10]
Figure 2005024354
Thereby, when the k-th eigenvalue is not degenerated, the k-th eigenvalue / eigenvector is obtained. If the eigenvalues are degenerated, using the Schmidt orthogonalization of the predetermined initial values with the already determined eigenvectors as the initial value of the inverse iteration method will reduce the parallelization effect, but the independent eigenvalues Eigenvectors can be obtained sequentially.
[0024]
As a third example, a method for suspending / resuming the parallel FDTD of the first example will be described with reference to FIGS. First, the collection integration processing unit 105 in FIG. If a signal or thread is used, an interruption request can be monitored without disturbing the calculation of FDTD. When the interruption request is detected, the process proceeds to the partial calculation synchronization process 320. Here, first, the collection integration processing unit 105 transfers an interruption request to each partial calculation processing unit 106. When the i-th partial calculation processing unit 106 detects this interruption request, the i-th partial calculation processing unit 106 advances the calculation of FDTD to a point where it is best and then pauses, and reports the processed simulation time ti to the collection integration processing unit 105. The collection integration processing unit 105 calculates the maximum value tmax of these simulation times and notifies each partial calculation processing unit 106. Each partial calculation processing unit 106 proceeds the FDTD processing as usual until the simulation time tmax. Here, the sound field collection integration processing 330 is performed. Each partial calculation processing unit 106 outputs the sound field (scalar field) p and medium velocity field (vector field) ν of the partial calculation region held at time tmax to the recovery integration processing unit 105. The collection integration processing unit 105 collects and integrates the data from each partial calculation processing unit 106, and creates the entire sound field p and medium velocity field ν at time tmax. Here, the sound field output processing 340 is performed. Here, the data of the sound field p, the medium velocity field ν, and the time tmax are stored in a predetermined file. This completes the interruption process. In the case of resumption, the initialization processing unit 101 may read the stored sound field p, medium velocity field ν, and time tmax data.
[0025]
【The invention's effect】
As described above, according to the present invention, a sound field can be calculated at high speed based on various acoustic models.
[Brief description of the drawings]
FIG. 1 is a configuration diagram of a sound field calculation apparatus using a parallel FDTD method.
FIG. 2 is an explanatory diagram of boundary data exchange between partial calculation processing units by a parallel FDTD method.
FIG. 3 is an explanatory diagram of interruption processing of a sound field calculation apparatus using a parallel FDTD method.
FIG. 4 is a configuration diagram of a sound field calculation apparatus using a parallel normal mode method.
FIG. 5 is a diagram showing a coordinate system describing an acoustic model.
[Explanation of symbols]
DESCRIPTION OF SYMBOLS 100 ... Structure of sound field calculation apparatus by parallel FDTD method and data flow, 101 ... Initialization processing unit, 102 ... Division processing unit, 103 ... Distribution processing unit, 104 ... Set of partial calculation processing unit, 105 ... Collection integration processing 106: Individual FDTD server in charge of partial calculation processing, 111 ... Data defining sound source, 112 ... Data defining medium field such as ocean (sound speed profile, depth profile, etc.), 113 ... Differentiation Lattice (range and step size, etc.), 114 ... Topology (region division method in partial calculation, mapping between partial calculation region and computer, etc.), 115 ... Difference scheme (coefficient of difference equation), 116 ... Initial sound field (and (Initial velocity field), 117 ... difference scheme divided for each partial calculation area, 118 ... initial sound field divided for each partial calculation area (O 119... Differenceing grid divided for each partial calculation area, 120... Startup data of partial calculation processing and collection integration processing, 121... Whole sound field, 200. ... partial calculation area 220 ... data exchange area 230 ... data exchange area thickness 300 ... sound field calculation interruption process 310 ... interrupt request detection process 320 ... partial calculation area synchronization process 330 ... sound field collection integration Processing: 340: Sound field output processing, 350: Sound field data at interruption, 400: Configuration and data flow of sound field calculation apparatus by parallel normal mode method, 401: Initialization processing unit, 402: Eigenvalue / eigenvector calculation processing 403... Collection integrated processing unit 411... Data defining a sound source, 412... Data defining a field of a medium such as the ocean (sound speed profile etc.), 41 ... (such as range and step size) grating for differencing, (coefficient of the difference of equations) 414 ... the difference scheme, 415 ... sound field.

Claims (5)

ＦＤＴＤ法の音場を並列計算する音場計算装置であって、
初期化処理部と分割処理部と配布処理部と部分計算処理部と統合処理部とからなり、
初期化処理部は媒質特性の場と格子化条件とを入力とし差分スキームおよび初期音場の格子データを出力し、
分割処理部は領域分割条件と上記差分スキームと初期音場とを入力とし
部分領域用に分割された差分スキームと初期音場と格子化条件とを出力し、
配布処理部は音源条件と領域分割条件と分割された差分スキームと初期音場と格子化条件とを入力とし
各部分計算処理および回収統合処理を起動して音源および上記入力データを該当するプロセスに出力し、
部分計算処理部はこれらのデータをもとに
所定の時間ステップごとに音場の更新および隣接する領域との必要な境界データの交換をおこない、
また音場を必要に応じて時間・空間方向にリサンプリングのうえ出力し、
統合処理部は各部分計算処理部からの出力を回収しつなぎ合わせてひとつの音場データを作成し・出力する
ことを特徴とする音場計算装置。A sound field calculation apparatus for calculating in parallel the sound field of the FDTD method,
It consists of an initialization processing unit, a split processing unit, a distribution processing unit, a partial calculation processing unit, and an integrated processing unit.
The initialization processing unit receives the medium characteristic field and the grid conditions, and outputs the difference scheme and the grid data of the initial sound field,
The division processing unit receives the region division condition, the above difference scheme, and the initial sound field, and outputs the difference scheme divided for the partial region, the initial sound field, and the lattice condition.
The distribution processing unit inputs the sound source condition, the region division condition, the divided difference scheme, the initial sound field, and the grid condition, and starts each partial calculation process and collection integration process, and converts the sound source and the input data into the corresponding process. Output,
Based on these data, the partial calculation processor updates the sound field and exchanges necessary boundary data with adjacent areas at predetermined time steps.
In addition, the sound field is output after resampling in the time and space directions as necessary.
A sound field calculation apparatus characterized in that the integrated processing unit collects and outputs the outputs from the respective partial calculation processing units to create and output one sound field data. 請求項１の音場計算装置において、
各部分計算処理部が境界データ交換用のバッファを備え、
複数の時間ステップごとにまとめて境界データの交換をおこなう
ことを特徴とする音場計算装置。The sound field calculation apparatus according to claim 1,
Each partial calculation processing section has a buffer for boundary data exchange,
A sound field calculation apparatus characterized in that boundary data is exchanged collectively at a plurality of time steps. 請求項１の音場計算装置において、
音場計算の中断要求を待ち受け、
中断要求があった場合各部分計算処理部を時刻をそろえたうえで止め、
その時点での各部分計算処理部の音場を回収・統合・保存して終了し、
また次回起動時に初期化処理部が上記で保存しておいた音場を取得し、
これにより音場計算の随時中断・再開を可能とする
ことを特徴とする音場計算装置。The sound field calculation apparatus according to claim 1,
Waiting for a request to interrupt the sound field calculation,
If there is an interruption request, stop each partial calculation processing unit after setting the time,
The sound field of each partial calculation processing unit at that time is collected, integrated, saved, and finished,
Also, when the next startup, the initialization processing unit acquires the sound field saved above,
This makes it possible to interrupt and resume sound field calculation at any time. ノーマルモード法の音場計算を並列的に実施する音場計算装置であって、
初期化処理部と複数の固有ベクトル算出部と回収統合部とからなり、
初期化処理部は音源、海洋、格子化の各条件を入力として
差分スキームを算出・出力し、
配布部は上記差分スキームを入力とし
そのコピーをそれぞれの固有ベクトル算出処理部に出力し、
固有ベクトル算出処理部は差分スキームを入力とし
これに対応する対称行列の所定の固有ベクトルを算出・出力し、
回収統合部は上記の固有ベクトルを入力とし
固有ベクトルを所定の方法で重ね合わせることにより音場を算出・出力する
ことを特徴とする音場計算装置。A sound field calculation device that performs sound field calculation of the normal mode method in parallel,
It consists of an initialization processing unit, multiple eigenvector calculation units, and a collection integration unit,
The initialization processing unit calculates and outputs the difference scheme with the sound source, ocean, and grid conditions as inputs.
The distribution unit takes the difference scheme as an input and outputs a copy thereof to each eigenvector calculation processing unit.
The eigenvector calculation processing unit receives a difference scheme as an input, calculates and outputs a predetermined eigenvector of a symmetric matrix corresponding to the difference scheme,
A sound field calculation apparatus characterized in that the recovery integration unit calculates and outputs a sound field by inputting the above eigenvectors and superimposing eigenvectors by a predetermined method. 請求項１〜４の音場計算装置において、
分散オブジェクトとの通信インタフェースを備え、
音源位置や観測位置などの条件を入力として随時受け付け、
その都度音場計算を実施または過去の音場計算結果を検索して
入力データに対応する音圧や伝播時間などの音場データを出力する
ことを特徴とする音場計算装置。In the sound field calculation apparatus according to claims 1 to 4,
It has a communication interface with distributed objects,
Accepts conditions such as sound source position and observation position as input
A sound field calculation apparatus that performs sound field calculation each time or searches past sound field calculation results and outputs sound field data such as sound pressure and propagation time corresponding to input data.
JP2003189156A 2003-07-01 2003-07-01 Parallel sound field calculator Pending JP2005024354A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2003189156A JP2005024354A (en) 2003-07-01 2003-07-01 Parallel sound field calculator

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2003189156A JP2005024354A (en) 2003-07-01 2003-07-01 Parallel sound field calculator

Publications (1)

Publication Number Publication Date
JP2005024354A true JP2005024354A (en) 2005-01-27

Family

ID=34187442

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2003189156A Pending JP2005024354A (en) 2003-07-01 2003-07-01 Parallel sound field calculator

Country Status (1)

Country Link
JP (1) JP2005024354A (en)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009077782A1 (en) * 2007-12-14 2009-06-25 University Of York Modelling wave propagation characteristics in an environment
WO2012101833A1 (en) * 2011-01-30 2012-08-02 富士通株式会社 Calculator and calculation method
CN110851781A (en) * 2019-11-01 2020-02-28 东南大学 Channel impulse response window estimation method based on frequency domain energy and group delay screening
CN115391069A (en) * 2022-10-27 2022-11-25 山东省计算中心（国家超级计算济南中心） Parallel communication method and system based on ocean mode ROMS
CN116244078A (en) * 2023-02-27 2023-06-09 青岛中海潮科技有限公司 Underwater sound field rapid calculation method based on multithreading and SIMD

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009077782A1 (en) * 2007-12-14 2009-06-25 University Of York Modelling wave propagation characteristics in an environment
WO2012101833A1 (en) * 2011-01-30 2012-08-02 富士通株式会社 Calculator and calculation method
CN110851781A (en) * 2019-11-01 2020-02-28 东南大学 Channel impulse response window estimation method based on frequency domain energy and group delay screening
CN110851781B (en) * 2019-11-01 2020-08-11 东南大学 Channel impulse response window estimation method based on frequency domain energy and group delay screening
CN115391069A (en) * 2022-10-27 2022-11-25 山东省计算中心（国家超级计算济南中心） Parallel communication method and system based on ocean mode ROMS
CN115391069B (en) * 2022-10-27 2023-02-03 山东省计算中心（国家超级计算济南中心） Parallel communication method and system based on ocean mode ROMS
CN116244078A (en) * 2023-02-27 2023-06-09 青岛中海潮科技有限公司 Underwater sound field rapid calculation method based on multithreading and SIMD
CN116244078B (en) * 2023-02-27 2023-12-01 青岛中海潮科技有限公司 Underwater sound field rapid calculation method based on multithreading and SIMD

Similar Documents

Publication Publication Date Title
Chen et al. Structural modal identification through high speed camera video: Motion magnification
US20180204341A1 (en) Ear Shape Analysis Method, Ear Shape Analysis Device, and Ear Shape Model Generation Method
WO2017031718A1 (en) Modeling method of deformation motions of elastic object
JP2018534712A (en) Unsupervised matching in fine-grained datasets for single view object reconstruction
EP4030346A1 (en) Information processing device and electronic apparatus equipped with same
CN112231847A (en) Transition position determination method and device, electronic equipment and storage medium
WO2017210228A1 (en) Digital camera with audio, visual and motion analysis
TW201812339A (en) Epicentral distance estimating device, epicentral distance estimating method, and computer-readable recording medium
CN112254798A (en) Method, system and medium for forecasting ocean vector sound field
Cai et al. Visualizing and quantifying small and nonstationary structural motions in video measurement
US20190178814A1 (en) State assessing device, state assessing method, and storage medium for storing program
JP2005024354A (en) Parallel sound field calculator
Forestier--Coste et al. Data-fitted second-order macroscopic production models
JP2015529814A (en) Methods for improving seismic horizon determination
Pires et al. Structural intensity assessment on shells via a finite element approximation
Li et al. Radial basis functions and level set method for image segmentation using partial differential equation
JP2010036600A (en) Device and method for estimating fluctuating pressure on hull surface by propeller and program
CN110136185B (en) Monocular depth estimation method and system
KR101892614B1 (en) Method and apparatus for aerodynamic analysis by manufacturing tolerance of blade
JPH0719818A (en) Three-dimensional movement predicting device
KR101436695B1 (en) Apparatus and method for detecting finger using depth image
Canelas et al. Topology design of inductors in electromagnetic casting using level-sets and second order topological derivatives
JP6685943B2 (en) Separation matrix design device, filter coefficient calculation device, method thereof, and program
US11317200B2 (en) Sound source separation system, sound source position estimation system, sound source separation method, and sound source separation program
JP2018064215A (en) Signal processing apparatus, signal processing method, and program