一种光路追踪计算方法及系统与流程

本发明涉及科学研究计算中的高性能计算领域，具体涉及一种光路追踪计算方法。
背景技术：
：激光在等离子体中的能量沉积问题是激光间接驱动惯性约束聚变(icf，inertialconfinementfusion)中的重要物理过程，其模拟精度会影响icf全过程数值模拟结果的正确性。在icf整体数值模拟中，目前激光主要采用几何光路建模，在模拟中需要通过求解光线与网格的相交关系来确定光线的运动轨迹。光路追踪算法有广泛的应用，例如无限通信系统、晶体光学、电磁辐射传播及碰撞检测、激光的传播等。在大规模并行机上实现快速高精度光路追踪计算一直是一件有挑战性的任务。计算光路与一系列(3-d或2-d)几何体的交点是光路追踪算法的核心组成部分。光的传播轨迹是由光线与几何体的交点组成的。相交点处及附近区域的物质特性决定了光线的状态和传播方向。如何提高算法的精度和效率往往因应用问题而异。在激光约束聚变问题中，激光作为聚变的源项出现在物理建模之中。具体地，通常采用几何光线近似来模拟激光在等离子体中吸收与传播。激光束被表示成多条激光光线(直线或曲线)。每条光线被赋予一定的激光能量。随着激光的传播，其携带的能量沉积到沿途的等离子体中。在每个计算网格的电子数密度线性分布的假设下，求交算法的核心就是求出光线参数方程与每个网格边线段的交点。在算法的设计上，一方面需要考虑到电子数密度为常数等可能的退化情形，另一方面需要考虑到计算机数字的有限表示位数，理论上相同数值上不同的算法可能会得到不同的数值解。由于浮点计算的局限性，很多常用的算法难以计算出理论上的精确解。而且在计算交点的算法中当判据或阈值达到计算机浮点精度的量级时，数值结果还可能会因计算机不同而不同，数值结果不同可能会引起光路传播方向的改变，或者传播路径长度的改变。ieee二进制浮点数算术标准(ieee754)是20世纪80年代以来最广泛使用的浮点数运算标准，为许多cpu与浮点运算器所采用，这个标准定义了表示浮点数的格式与反常值，一些特殊数值，例如：无穷(inf)与非数值(nan)，以及这些数值的“浮点数运算符”；它也指明了四种数值舍入规则和五种例外状况，包括例外发生的时机与处理方式。因为ieee754标准指明了浮点数的有效存储位数，导致浮点数在计算机中运算结果的不精确性来自于浮点数的有限位数无法精确表达自然界中的一个自然数，进一步导致了一些算法在理论上正确无误，但在计算机上则可能计算结果与真实数值不同，而计算错误。技术实现要素：为解决上述技术的不足，本发明提供一种光路追踪计算方法，根据浮点稳定性计算理论，采用扩展精度浮点数计算方法求解光线参数方程，解决了减法相消带来的计算错误，提高了光线参数方程求解的准确性，能精准地实现光路的追踪。一种光路追踪计算方法，其改进之处在于，所述方法包括：采集被测激光的参数值；基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；基于所述激光传播时间计算所述激光的轨迹点；所述光线参数方程包括一元四次方程。优选地，所述基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：将一元四次光线参数方程简化为一元三次方程；基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；其中，所述扩展精度浮点数包括精度高于双精度的浮点数。优选地，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程包括：计算所述一元三次方程的判别式；基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。优选地，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。优选地，所述基于所述求解结果确定所述一元三次方程的根包括：采用牛顿迭代法对所述求解结果进行迭代计算；当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；将所述迭代计算结果作为所述一元三次方程的根。优选地，所述是否出现减法相消按下式进行判断：式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。优选地，所述一元三次方程的判别式如下式所示：式中：δ：判别式；p、q：消去二次项后的一元三次方程的系数；其中：p按下式进行计算：q按下式进行计算：优选地，所述一元三次方程的判别式＜0时，所述一元三次方程的根的个数为3个，求根公式如下式所示：式中：x1、x2、x3：一元三次方程的根；θ：采用三角函数法求解一元三次方程时的参数；其中，θ按下式进行计算：优选地，所述基于所述激光传播时间计算激光的轨迹点包括：将所述激光传播时间代入到关于光线传播时间的光线参数方程，得到激光的轨迹点。一种光路追踪计算系统，包括参数采集模块、时间获取模块和轨迹点获取模块；参数采集模块：用于采集被测激光的参数值；时间获取模块：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；轨迹点获取模块：基于所述激光传播时间计算所述激光的轨迹点；所述光线参数方程包括一元四次方程。与最接近的已有技术比，本发明提供的技术方案具有以下有益效果：本发明提供的技术方案，根据浮点稳定性计算理论，采用扩展精度浮点数计算方法求解光线参数方程，解决了减法相消带来的计算错误，提高了光线参数方程求解的准确性，精准地实现光路的追踪。本发明提供的技术方案，采用扩展精度运算、迭代提炼精度和重根检测的方法优化解的精度，提高了光线参数方程求解算法的稳定性。附图说明图1是本发明光路追踪计算方法的示意图；图2是本发明基于扩展精度浮点数采用三角函数法求解一元三次方程的流程图；图3是本发明光路追踪计算系统的示意图。具体实施方式为了更好地理解本发明，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。实施例一、一种光路追踪计算方法，如图1所示，包括：步骤1：采集被测激光的参数值；步骤2：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；步骤3：基于所述激光传播时间计算所述激光的轨迹点；所述光线参数方程包括一元四次方程。步骤1：采集被测激光的参数值。步骤2：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间。具体地，所述采用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：将一元四次光线参数方程简化为一元三次方程；基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；其中，所述扩展精度浮点数包括精度高于双精度的浮点数。具体地，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程包括：计算所述一元三次方程的判别式；基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。具体地，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。具体地，所述基于所述求解结果确定所述一元三次方程的根包括：采用牛顿迭代法对所述求解结果进行迭代计算；当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；将所述迭代计算结果作为所述一元三次方程的根。具体地，所述减法相消按下式进行判断：式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。具体地，所述一元三次方程的判别式如下式所示：式中：δ：判别式；p、q：消去二次项后的一元三次方程的系数；其中：p按下式进行计算：q按下式进行计算：具体地，所述一元三次方程的判别式＜0时，所述一元三次方程的根的个数为3个，求根公式如下式所示：式中：x1、x2、x3：一元三次方程的根；θ：采用三角函数法求解一元三次方程时的参数；其中，θ按下式进行计算：步骤3：所述基于所述激光传播时间计算激光的轨迹点。将所述激光传播时间代入到关于光线传播时间的光线参数方程，得到激光的轨迹点。实施例二、一种光路追踪计算方法，如图1所示，所述方法包括：步骤1：采集被测激光的参数值；步骤2：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；步骤3：基于所述激光传播时间计算所述激光的轨迹点；所述光线参数方程包括一元四次方程。步骤1：采集被测激光的参数值。在几何光学近似下，光线传播方程可以描述为其中是光线的位置矢量，ω和分别是激光的频率和传播速度t表示激光传播的时间，c是真空中的光速，ne是电子数密度，临界电子数密度nc＝ω2me/4πe2，me是电子质量，e为电子电量。描述激光在等离子体中传播的是麦克斯韦方程，它的求解要求网格尺度是激光波长量级的，在时间尺度方面也有同样的要求。而目前的模拟程序还远远不能达到。在这种情况下提出了几何光学近似下的光路追踪法。它的基本思想是：将激光束分成若干条光线，每条光线初始时赋予一定的激光能量，在等离子体中传播，其行为满足光线参数方程，沿着光线轨迹，激光能量通过逆韧致过程和其他过程被等离子体吸收。由于激光传播的速度很快，可以假设在一个流体力学步长内能量沉积是瞬时的，在一条光线传播过程中，等离子体的状态是不变的。为将三维光线应用于以z为对称轴的二维柱几何数值计算中，即其中(kr，kθ，kz)表示光线在柱坐标系下的速度，可通过坐标变换将直角坐标系中的光线参数方程变换为柱坐标系下的光线参数方程：通过无量纲化略去关于网格的附标，并且假定在每一个网格上电子数密度的变化是线性的，则可以得到柱坐标系中每个二维网格上的三维光线参数方程：其中，(r0，z0，θ0)与分别为光线进入点的位置与波矢。求解上述方程可以得到关于t的光线参数方程，具体地，a)当h*＝0时，光线参数方程是：b)当gr＝0时，光线参数方程是：c)当h*≠0且gr≠0时，无法精确求解，记令s＝y′(r)，则光线参数方程可近似表示为：基于以上光线参数方程的具体形式，可以计算光线参数方程与二维柱几何网格的交点。具体地，要分别计算光线与每条网格边的交点。把网格边所在的直线方程记作r＝az+b，a＝(r2-r1)/(z2-z1)，b＝r1-az1，其中(r1，z1)和(r2，z2)分别为任一网格边的两个端点，由此可以相应地得到如下关于t的一元二次方程和一元四次方程。a)当h*＝0时，b)当gr＝0时，c)当h*≠0且gr≠0时，求解出上述一元二次方程或一元四次方程就可以得到参数t的值，代入光线参数方程即可得到交点的坐标值。求解一元二次或四次方程可能会得到多个根，判断真解的条件是t＞10-6且交点在网格边上。如有多个根满足真解条件，应选取最小值作为光线的轨迹点。根据光线的轨迹点就可以求出光线穿过网格的路程，继而计算出这条光线沉积的激光能量。费拉里算法是求解一元四次方程的一种常用算法，具体如下：对于一元四次方程a′x4+b′x3+c′x2+d′x+e＝0将方程的最高阶次项化为1，执行将方程移项，化为x4+b′x3＝-c′x2-d′x-e方程两边同时加上使左边配成完全平方式，如下所示：引入参数y并解出y，使得方程左右两边都配成完全平方式。首先，左右都加上得到需要将右边配成完全平方式，即右边的判别式为0，如下式所示：化简为一个三次方程-y3+c′y2+(4e-b′d′)y+d′2+eb′2-4ec′＝0然后求解这个一元三次方程，会得到y的三个值，取其中任何一个即可。可以得到原方程等价于可以得到两个一元二次方程对这两个一元二次方程化简求解这两个一元二次方程，即可得到一元四次方程的四个实根。在上述费拉里算法中，需要求解一个一元三次方程，该方程的解是否精确将直接影响费拉里算法解出的一元四次方程根的正确性。步骤2：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间。具体地，所述采用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：将一元四次光线参数方程简化为一元三次方程；基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；其中，所述扩展精度浮点数包括精度高于双精度的浮点数。计算机中常用的浮点数是单精度，双精度浮点数，本发明使用了扩展精度浮点数，这些类型浮点数的指数和有效数位数如下：表1浮点数的位数分布类型占用位数指数位数小数位数单精度32823双精度641152如图2所示，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程包括：计算所述一元三次方程的判别式；基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。具体地，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。具体地，所述基于所述求解结果确定所述一元三次方程的根包括：采用牛顿迭代法对所述求解结果进行迭代计算；当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；将所述迭代计算结果作为所述一元三次方程的根。牛顿迭代法求解方程f(x)＝0。给定初始估计解x0，牛顿迭代法是通过以下公式迭代求得下一步解，如下式所示：牛顿迭代法能够收敛到解的精度的极限取决于能计算到的精确程度。在浮点计算中，减法相消发生在两个数值上非常接近的浮点数相减，这经常会导致大量的有效数字丢失。考虑一个函数其值域为[0，0.5]，在x＝1.2×10-5，有效数字为10位的情况下cosx＝0.99999999991-cosx＝0.0000000001这明显计算错误，这个问题的原因是1-cosx的有效位数只有一位，1-cosx的浮点计算结果是精确的，但减法操作把浮点误差放大了。在这个函数中，利用恒等式cosx＝1-2sin2(x/2)，将方程重写为以下形式从而避免舍入误差假设有两个相近的浮点数m和n，有和和是m和n欲代表自然数，δm和δn是数据中的舍入误差，或者之前计算的误差，计算x＝m-n，有当|m-n|＜＜|m|+|n|时，结果和自然界中对应的结果会产生严重的相对误差。但是，当δm＝δn＝0，即浮点数能精确表示自然数时，不会产生减法相消误差。所述减法相消按下式进行判断：式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。所述一元三次方程如下式所示：ax3+bx2+cx+d＝0。具体地，所述一元三次方程的根的个数为3个时，求根公式如下式所示：式中，p、q：消去二次项后的一元三次方程的系数；此时，判别式小于零。其中，采用三角函数法求解一元三次方程时的参数θ按下式进行计算：p按下式进行计算：判别式按下式进行计算：其中：δ：判别式；q按下式进行计算：若判别式等于零：若p＝q＝0，一元三次方程有一个根，求根公式如下式所示：否则一元三次方程有两个根：若q＜0，求根公式如下式所示：否则求根公式如下式所示：若判别式＞0，一元三次方程有一个根：若求根公式如下式所示：否则求根公式如下式所示：步骤3：将所述激光传播时间代入光线参数方程，得到光线的轨迹点包括：将所述激光传播时间代入到关于光线传播时间的光线参数方程，得到激光的轨迹点。表2列出了在光路追踪高精度计算方法中使用判别式扩展精度、求根公式减法相消检测与扩展精度、牛顿迭代法、重根检测与重启牛顿迭代的各种情况，以及各种情况下经过对1300万个一元三次方程进行计算，统计的求解正确率。由此可见使用本发明提供的技术方案可获得100％正确率。而经典三角函数算法，没有使用判别式扩展精度、求根公式减法相消检测与扩展精度、牛顿迭代法、重根检测与重启牛顿迭代的情况下对这些方程的正确率约为85％。表2比较表实施例三、一种光路追踪计算系统，如图3所示，包括参数采集模块、时间获取模块和轨迹点获取模块；参数采集模块：用于采集被测激光的参数值；时间获取模块：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；轨迹点获取模块：基于所述激光传播时间计算所述激光的轨迹点；所述光线参数方程包括一元四次方程；具体地，轨迹点获取模块中基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：将一元四次光线参数方程简化为一元三次方程；基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；其中，所述扩展精度浮点数包括精度高于双精度的浮点数。具体地，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程：计算所述一元三次方程的判别式；基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。具体地，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。具体地，所述基于所述求解结果确定所述一元三次方程的根包括：采用牛顿迭代法对所述求解结果进行迭代计算；当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；将所述迭代计算结果作为所述一元三次方程的根。具体地，所述减法相消按下式进行判断：式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。具体地，所述一元三次方程的判别式如下式所示：式中：δ：判别式；p、q：消去二次项后的一元三次方程的系数；其中：p按下式进行计算：q按下式进行计算：具体地，所述一元三次方程的判别式＜0时，所述一元三次方程的根的个数为3个，求根公式如下式所示：式中：x1、x2、x3：一元三次方程的根；θ：采用三角函数法求解一元三次方程时的参数；其中，采用三角函数法求解一元三次方程时的参数θ按下式进行计算：本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员依然可以对本发明的具体实施方式进行修改或者等同替换，这些未脱离本发明精神和范围的任何修改或者等同替换，均在申请待批的本发明的权利要求保护范围之内。当前第1页12