CN111967148A - 一种粒子模拟的Voronoi图粒子合并算法 - Google Patents
一种粒子模拟的Voronoi图粒子合并算法 Download PDFInfo
- Publication number
- CN111967148A CN111967148A CN202010757434.1A CN202010757434A CN111967148A CN 111967148 A CN111967148 A CN 111967148A CN 202010757434 A CN202010757434 A CN 202010757434A CN 111967148 A CN111967148 A CN 111967148A
- Authority
- CN
- China
- Prior art keywords
- particle
- particles
- merging
- voronoi
- simulation
- 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.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及粒子模拟数值模拟领域,具体为一种基于粒子模拟的Voronoi图粒子合并算法。本发明在合并前对粒子列表进行预处理,筛掉权重较大的粒子,避免因合并后粒子权重过大而引起较大的数值噪声,同时,提升了整个计算过程的速度。并通过预处理规则的选取,使得期望权重的选取更佳合理,从而达到快速找出相似粒子分组的目的,提高了Voronoi图粒子合并算法的计算效率。本发明克服了PIC方法模拟过程中粒子数目过大的限制,该算法在降低了模拟粒子规模同时保证了合并前后物理特性的不变,极大地提高了计算效率,并且可以方便的调节用户输入参数来控制最终合并误差和用户期望的剩余粒子数。
Description
技术领域
本发明涉及粒子模拟数值模拟领域,具体为一种基于粒子模拟的Voronoi图粒子合并算法,提高粒子模拟的数值计算效率。
背景技术
粒子模拟(PIC)方法是通过跟踪大量在自洽和外加电磁场作用下的带电粒子的运动来研究粒子宏观特性及运动规律的动力学模拟方法。PIC方法的基本流程为:根据粒子的位置与速度计算得到网格节点上的电荷密度和电流密度;将此电荷和电流密度代入泊松方程或麦克斯韦方程计算出网格节点上的电场与磁场;根据电场与磁场求得粒子所在位置上的洛伦兹力,代入粒子运动方程得到新时刻的粒子位置与速度。重复上述流程直至模拟时间结束或达到收敛条件。
PIC方法通过跟踪大量带电粒子的运动,根据感兴趣的问题对某些物理量做统计平均,可以较为准确的描述粒子的宏观物理特性和微观运动过程。目前该方法已经在等离子体的许多领域得到了广泛的应用,例如空间物理研究、热核聚变、惯性约束聚变等。但是PIC方法模拟物理问题时受稳定性条件约束的影响,时间步长往往取值较小,而许多物理问题的模拟时间相较于PIC的时间步长跨度较大,并且PIC每个时间步都要跟踪大量的带电粒子运动以及求解粒子与电磁场间的相互作用,因此计算量十分巨大。同时,采用PIC方法模拟诸如电离、二次电子倍增等带电粒子数目不断增加的物理问题时,不可避免地会产生庞大的数值计算负担,从而导致模拟难以为继。
为了克服粒子数目过多对PIC模拟造成的限制,需要引入粒子合并方法,其中Voronoi图的粒子合并算法在保证模拟问题物理特性不变的同时可以有效的减少粒子数目,进而提高模拟效率。Voronoi图的粒子合并算法主要分为三个步骤,首先是初始化用户参数;其次是构造Voronoi图数据结构对待合并的粒子进行分组;最后是对分组后的粒子进行合并。通过构造Voronoi图可以高效的找到属性“接近”的粒子,对属性接近的粒子进行合并,可以保证合并前后粒子的位置和速度分布特性不会发生明显的变化。通过实施不同的合并方案,如动量守恒或能量守恒等方式对分组完成的粒子进行合并,可以适用于不同的物理问题。但是传统的Voronoi图粒子合并算法未能考虑粒子的权重变化差异。当待合并的粒子权重变化较大时,权重较高的粒子参与合并,会造成合并后的部分粒子权重过大,进而引起较大噪声,严重的影响最终模拟结果的准确性。同时粒子权重变化较大,会使得Voronoi图的构造难以收敛,不能划分出合适的粒子分组,从而降低Voronoi图粒子合并算法的计算效率。
发明内容
针对粒子模拟中现有Voronoi图的粒子合并技术存在的问题或不足,本发明提供了一种粒子模拟的Voronoi图粒子合并算法。该方法通过引入预处理技术,筛掉权重较高的粒子,找到合适的待合并粒子进行合并,从而减少了合并后的数值噪声和提高计算效率。
具体技术方案如下:
步骤1、初始化用户参数。设置用户参数T、ratio和Nmin,其中T是最大合并粒子数目,即不超过T个粒子可以合并成一个粒子,避免因过多粒子参与融合造成合并后的粒子权重过大的情况,ratio是剩余粒子期望比率,用于决定用户期望的剩余粒子数,Nmin是粒子合并阈值,用于判断当前网格是否需要进行粒子合并。
步骤2、预处理粒子列表,筛选出待合并粒子列表。
设当前粒子列表数目为N。若N>Nmin,则需要进行粒子合并,否则不需要进行合并。通过ratio计算用户期望剩余粒子数Nopt=ratio*N,进而得到期望权重Wopt=max(1,N/Nopt)。
将所有粒子的权重与期望权重进行比较:对于权重为Wi的第i个粒子,若Wi<Wopt,则将该粒子加入待合并粒子列表;若Wi≥Wopt,则该粒子不参与合并,直接保留其原有属性。以期望权重为基准,筛选出权重较低的粒子进行合并,避免了权重较高的粒子参与合并而引起较大噪声。
进一步的,所述步骤2中的合并,以合并后粒子的权重不超过期望权重为基准,合并后的粒子循环参与合并。现有技术中使用固定的期望权重,在不同的求解时刻,不能筛选出合适的参与合并的粒子。而本发明通过参数ratio求得的期望权重,每个时刻都可以得到一个更加合适的期望权重作为筛选基准,降低了模拟噪声从而提高计算的准确性;并且合并后粒子会再参与合并,进行多次合并后会更进一步的提升计算效率。
步骤3、通过构造Voronoi图对待合并的粒子进行分组;
构造Voronoi图的过程:将模拟区域划分成一个个的Voronoi单元(简称V单元),属性相似的粒子被分配到各自对应的V单元内,完成待合并粒子的分组。
进一步的,具体分组步骤如下:
步骤3.2、将和分别与容差系数errx和errp进行比较,其中errx和errp表示V单元内粒子位置和速度的偏离程度判定系数,根据模拟误差的需要(不超过模拟误差)进行取值。errx和errp值越小,表明粒子越“接近”,反之表明粒子较为离散。若所有方向都有和则表明当前V单元内的粒子足够接近(满足合并条件),停止分裂并标记为true。
对于需要分裂的V单元,找到位置和速度偏离最大的维度方向k。将V0按照垂直于k方向进行分裂,则V0分裂成两个新的V单元,记为V1和V2。将V0内的粒子信息分类装入V1和V2(例如在位置的k方向进行划分,则将xi≤X0,k的粒子放入V1,xi>X0,k的粒子放入V2),分类完成后删除V0。
步骤3.3、对V1和V2执行步骤3.1和3.2,直至所有的V单元都遍历完成。此时Voronoi图构造完成,粒子被划分到各个V单元中,且每个V单元内的粒子属性最为接近。
步骤4、实施不同的合并方案对分组粒子进行合并。对于划分好的各个V单元,标记为false的V单元内的粒子直接装入合并后粒子列表。标记为true的V单元(包含Vnum个粒子),采用已有的能量守恒或动量守恒的方式按照(Vnum:1)进行合并。所有V单元合并完成后,当前网格合并结束,遍历模拟区域的所有网格直至合并完成。
本发明与现有用于粒子模拟的粒子合并方法相比具有以下优势,该算法可以方便的调节用户输入参数来控制最终合并误差和用户期望的剩余粒子数。在合并前对粒子列表进行预处理,筛掉权重较大(对于大于期望权重的认定为权重过大)的粒子,仅满足需求的粒子才会参与合并,避免因合并后粒子权重过大而引起较大的数值噪声,同时,提升了整个计算过程的速度。并进一步的通过预处理规则的选取,使得期望权重的选取更佳合理,从而提升了V图的构造速度(期望权重过大会导致V图构造时生成过多的V单元,从而构造过程变慢),从而达到快速找出相似粒子分组的目的,提高了Voronoi图粒子合并算法的计算效率。
综上所述,本发明可以克服PIC方法模拟过程中粒子数目过大的限制,该算法在降低了模拟粒子规模同时保证了合并前后物理特性的不变,极大地提高了计算效率。
附图说明
图1为二维轴对称放电室模型计算实例示意图;
图2为合并前后粒子的轨迹分布对比图;
图3为合并前后粒子各方向速度分量分布对比图;
图4为合并前后粒子的数目变化对比图。
具体实施方式
下面以一个离子推进器放电室的磁约束为例,通过具体实施过程对本发明作进一步详细说明。
本实施例模拟了一个二维轴对称离子推进器放电室模型,其外围包裹着的三个钐钴永磁体环,分别放置于放电室的前壁面、侧壁面和后壁面,结构示意图如图1所示。
采用PIC法模拟粒子的运动,观测粒子在磁场约束下的运动轨迹。放电室半径为20cm,长18cm,空间步长Δr=Δz=2mm,时间步长Δt=2×10-12s。初始时刻将粒子随机分布在4cm<r,z<14cm区域内,每个网格包含100个粒子,每个粒子的权重在[1,50]中随机分配,粒子的速度采用Maxwell分布。模拟运行时长60ns,合并周期为2ns。
步骤1、初始化用户参数T=5、ratio=0.5和Nmin=15。因此放电室内的粒子最多可以进行(5:1)的合并,避免了合并后的粒子权重过大。通过ratio计算得到用户期望的剩余粒子数。设定Nmin=15,符合PIC模拟放电室的网格平均粒子数目要求。
步骤2、随着模拟的进行,当模拟时间达到合并周期时启动粒子合并程序。假设当前网格包含的粒子数目为N。比较N与Nmin,判断是否需要进行粒子合并,若不需要合并则跳过该网格,否则计算用户期望剩余粒子数Nopt=ratio*N,求解期望权重Wopt=max(1,N/Nopt)。遍历网格内的粒子列表,对于Wi<Wopt的粒子加入到待合并粒子列表中,对于Wi>=Wopt的粒子,权重过大,不参与合并,直接保留。只要粒子的权重不超过期望权重,合并后的粒子就将循环参与合并。
步骤3、通过待合并粒子列表构建Voronoi图,对属性相似的粒子进行分组。
步骤3.2、将和分别与容差系数errx和errp进行比较来决定是否停止划分,当前的位置容差errx=4e-2,动量容差errp=4e-6。较小的errx和errp会划分出较多的V单元,消耗更多的时间,但同时找到的粒子也越接近,合并后的粒子速度分布也与原始系统更为接近。若所有方向都有和则表明当前V单元内的粒子足够接近(满足合并条件),停止分裂并标记为true。
若有任意方向的或不满足条件,则比较当前V单元内的粒子数目是否小于T,若小于T,则表明V单元足够小,此时V单元停止分裂标记为false,否则当前V单元需要继续分裂。对于需要分裂的V单元,找到位置和速度偏离最大的维度方向k。将V0按照垂直于k方向进行分裂,则V0分裂成两个新的V单元,记为V1和V2。将V0内的粒子信息分类装入V1和V2(例如在位置的k方向进行划分,则将xi≤X0,k的粒子放入V1,xi>X0,k的粒子放入V2),分类完成后删除V0。
步骤3.3、对V1和V2执行步骤3.1和3.2,直至所有的V单元都遍历完成,至此Voronoi图构建完成。Voronoi图的构建根据粒子的位置和动量信息将网格划分成多个V单元,每个V单元内的粒子即为属性相似的粒子。
步骤4、划分好V单元后,对于标记为false的V单元内的粒子直接保留,装入到合并后的粒子列表,对于被标记为true的V单元,采用能量守恒或动量守恒的方式按照(Vnum:1)进行合并。此处采用能量守恒的方式进行合并。假设当前V单元内权重为w1,w2…wn,位置为x1,x2…xn以及速度为υ1,υ2…υn的粒子需要合并,则合并后的粒子新位置xnew为
粒子加权平均速度υavg为
合并后的粒子新速度υnew为
比较有无粒子合并情况下,60ns时粒子的相空间分布情况。图2为粒子的轨迹分布对比图,图3为粒子的三个方向速度分量的相空间分布对比图,图4为合并前后粒子的数目变化对比图。通过对比发现,使用本发明的粒子合并方法,在进行多次合并后,粒子的位置和速度相空间分布基本保持不变,粒子规模有较大程度的减小,在提高模拟效率的同时保证了模拟问题的物理特性一致。
Claims (4)
1.一种粒子模拟的Voronoi图粒子合并算法,其特征在于,包括以下步骤:
步骤1、初始化用户参数:设置用户参数T、ratio和Nmin,其中T是最大合并粒子数目,即不超过T个粒子可以合并成一个粒子,ratio是剩余粒子期望比率,用于决定用户期望的剩余粒子数,Nmin是粒子合并阈值,用于判断当前网格是否需要进行粒子合并。
步骤2、预处理粒子列表,筛选出待合并粒子列表:
设当前粒子列表数目为N,若N>Nmin,则需要进行粒子合并,否则不需要进行合并,通过ratio计算用户期望剩余粒子数Nopt=ratio*N,进而得到期望权重Wopt=max(1,N/Nopt);
将所有粒子的权重与期望权重进行比较:对于权重为Wi的第i个粒子,若Wi<Wopt,则将该粒子加入待合并粒子列表;若Wi≥Wopt,则该粒子不参与合并,直接保留其原有属性。
步骤3、通过构造Voronoi图对待合并的粒子进行分组:
将模拟区域划分成一个个的Voronoi单元,属性相似的粒子被分配到各自对应的V单元内,完成待合并粒子的分组。
步骤4、实施合并方案对分组粒子进行合并:
对于划分好的各个Voronoi单元,标记为false的Voronoi单元内的粒子直接装入合并后粒子列表;标记为true的Voronoi单元,按合并方案进行合并,所有Voronoi单元合并完成后,当前网格合并结束,遍历模拟区域的所有网格。
2.如权利要求1所述粒子模拟的Voronoi图粒子合并算法,其特征在于:所述步骤4中的合并方案为能量守恒或动量守恒的合并方案。
3.如权利要求1所述粒子模拟的Voronoi图粒子合并算法,其特征在于:所述步骤2中的合并,以合并后粒子的权重不超过期望权重为基准,合并后的粒子循环参与合并。
4.如权利要求1所述粒子模拟的Voronoi图粒子合并算法,其特征在于,所述步骤3中构造Voronoi图的具体步骤如下:
步骤3.2、将和分别与容差系数errx和errp进行比较,其中errx和errp表示Voronoi单元内粒子位置和速度的偏离程度判定系数,不超过模拟误差;若所有方向都有和则表明当前Voronoi单元内的粒子满足合并条件,停止分裂并标记为true;
若有任意方向的或不满足条件,则比较当前Voronoi单元内的粒子数目是否小于T,若小于T,则表明Voronoi单元足够小,此时Voronoi单元停止分裂标记为false;否则当前Voronoi单元需要继续分裂;
对于需要分裂的Voronoi单元,找到位置和速度偏离最大的维度方向k;将V0按照垂直于k方向进行分裂,则V0分裂成两个新的Voronoi单元,记为V1和V2;将V0内的粒子信息分类装入V1和V2,分类完成后删除V0;
步骤3.3、对V1和V2执行步骤3.1和3.2,直至所有的Voronoi单元都遍历完成。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010757434.1A CN111967148B (zh) | 2020-07-31 | 2020-07-31 | 一种粒子模拟的Voronoi图粒子合并算法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010757434.1A CN111967148B (zh) | 2020-07-31 | 2020-07-31 | 一种粒子模拟的Voronoi图粒子合并算法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111967148A true CN111967148A (zh) | 2020-11-20 |
CN111967148B CN111967148B (zh) | 2023-07-07 |
Family
ID=73363729
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010757434.1A Active CN111967148B (zh) | 2020-07-31 | 2020-07-31 | 一种粒子模拟的Voronoi图粒子合并算法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111967148B (zh) |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110219344A1 (en) * | 2010-03-05 | 2011-09-08 | International Business Machines Corporation | Spatial Correlation-Based Estimation of Yield of Integrated Circuits |
CN103440163A (zh) * | 2013-09-09 | 2013-12-11 | 中国科学院近代物理研究所 | 使用gpu并行实现的基于pic模型的加速器仿真方法 |
CN106855900A (zh) * | 2015-12-08 | 2017-06-16 | 三星电子株式会社 | 基于流体粒子建模空气泡的运动的方法和设备 |
CN107705321A (zh) * | 2016-08-05 | 2018-02-16 | 南京理工大学 | 基于嵌入式系统的运动目标检测与跟踪方法 |
CN108416107A (zh) * | 2018-02-05 | 2018-08-17 | 电子科技大学 | 一种应用于pic的推动粒子运动有限元算法 |
CN109063332A (zh) * | 2018-08-02 | 2018-12-21 | 中国科学院自动化研究所 | 重心Voronoi图的规整性提升方法 |
CN110309248A (zh) * | 2019-06-26 | 2019-10-08 | 东南大学 | 一种基于Voronoi图的交通道路网络自动划分交通小区的方法 |
-
2020
- 2020-07-31 CN CN202010757434.1A patent/CN111967148B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110219344A1 (en) * | 2010-03-05 | 2011-09-08 | International Business Machines Corporation | Spatial Correlation-Based Estimation of Yield of Integrated Circuits |
CN103440163A (zh) * | 2013-09-09 | 2013-12-11 | 中国科学院近代物理研究所 | 使用gpu并行实现的基于pic模型的加速器仿真方法 |
CN106855900A (zh) * | 2015-12-08 | 2017-06-16 | 三星电子株式会社 | 基于流体粒子建模空气泡的运动的方法和设备 |
CN107705321A (zh) * | 2016-08-05 | 2018-02-16 | 南京理工大学 | 基于嵌入式系统的运动目标检测与跟踪方法 |
CN108416107A (zh) * | 2018-02-05 | 2018-08-17 | 电子科技大学 | 一种应用于pic的推动粒子运动有限元算法 |
CN109063332A (zh) * | 2018-08-02 | 2018-12-21 | 中国科学院自动化研究所 | 重心Voronoi图的规整性提升方法 |
CN110309248A (zh) * | 2019-06-26 | 2019-10-08 | 东南大学 | 一种基于Voronoi图的交通道路网络自动划分交通小区的方法 |
Non-Patent Citations (1)
Title |
---|
郭胜龙: "离子推进器放电室电磁特性的研究", 中国博士学位论文全文数据库 基础科学辑 * |
Also Published As
Publication number | Publication date |
---|---|
CN111967148B (zh) | 2023-07-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Kaveh et al. | A novel heuristic optimization method: charged system search | |
CN107067136B (zh) | 电动汽车充电分配方法及装置 | |
CN104732044B (zh) | 基于差分进化算法的多层频率选择表面复合吸波结构及材料的优化设计方法 | |
Pei et al. | Triple and quadruple comparison-based interactive differential evolution and differential evolution | |
Ohnishi et al. | Correlation between potential well structure and neutron production in inertial electrostatic confinement fusion | |
O'Brien et al. | Automatic simplification of particle system dynamics | |
Tang et al. | A new Nash optimization method based on alternate elitist information exchange for multi-objective aerodynamic shape design | |
Procassini et al. | Particle simulation model of transport in a bounded, Coulomb collisional plasma | |
Nakayama et al. | Numerical simulation of ion beam optics for multiple-grid systems | |
CN110766125A (zh) | 一种基于人工鱼群算法的多目标武器——目标分配方法 | |
CN111967148A (zh) | 一种粒子模拟的Voronoi图粒子合并算法 | |
CN111709183B (zh) | 一种基于遗传算法的中子管加速系统优化方法 | |
Park | Guiding flows for controlling crowds | |
CN116485043B (zh) | 一种翼伞集群系统归航多目标优化方法 | |
Yu et al. | A QPSO algorithm based on hierarchical weight and its application in cloud computing task scheduling | |
Muth et al. | Contacts between many bodies | |
CN113392568B (zh) | 飞行器气动特性并行模拟中动态计算域的负载平衡方法 | |
CN114841079A (zh) | 一种基于深度神经网络优化的连续碰撞检测方法 | |
CN115618699A (zh) | 一种适用于超声速稀薄流动模拟的优化型信息保存法 | |
CN113987898B (zh) | 一种放电室中气体放电的粒子模拟加速方法 | |
CN114123216A (zh) | 新能源电力系统的经济负荷分配方法、系统及介质 | |
Boronina et al. | The parallel threedimensional PIC code for the numerical modeling of ultrarelativistic beams | |
Blanco et al. | Paul’s Trap | |
Meyer et al. | Asynchronous cycling as a convergence acceleration method in particle simulation of direct current glow discharges | |
Andersson et al. | Response surface methods and Pareto optimization in crashworthiness design |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |