精华内容
下载资源
问答
  • 建立二维、三维单胞模型的周期性边界条件python脚本程序
  • 建立二维、三维单胞模型的周期性边界条件python脚本程序
  • 周期性边界条件

    2021-04-06 16:17:52
    其三就是周期性边界条件, 一定要想明白, 不管是写参数还是写程序都得经常敲敲脑瓜子. 今天我们就来谈谈这个periodic boundary condition. 我们直到今天, 仍然只能模拟一个很小的体系. 一方面, 哪怕是前两天报道了一...

    在我刚开始做模拟的时候, 我老板说别的你可以不知道, 有几个原则必须烙在你的脑子里, 稍一不注意就会完犊子. 其一是我们能模拟的体系\元胞\盒子非常非常小, 比你想象中的最小的还要小; 其二是我们能模拟的时间跨度非常非常短, 分析数据的时候一定要提醒自己; 其三就是周期性边界条件, 一定要想明白, 不管是写参数还是写程序都得经常敲敲脑瓜子.

    今天我们就来谈谈这个periodic boundary condition.

    我们直到今天, 仍然只能模拟一个很小的体系. 一方面, 哪怕是前两天报道了一个对有着两亿个原子的动力学模拟, 实际上还是没有质的变化; 另一方面即便我们可以倾尽家产, 包下太湖或者微光, 可以, 但没必要. 我们通常认为, 对于一个极小元胞的模拟, 就能反映出整体的特征趋势, 剩下的工作应该用你聪明的大脑来展开.

    考虑一个极小的立方体液滴单元, 里面有一千个原子. 其中有 个在内部, 可以看作完全处于液体环境中, 而剩下的接近一半, 都在这个单元的表皮上, 显然不会有着相同的性质. 即使是到了 个原子这样一个量级, 仍然有接近6%有着不同的性质. 而我们所取得这个液体单元, 必须要是对宏观普适的抽样, 那么怎么做呢?

    这个问题可以通过由Born and von Karman 在1912年提出的周期性边界条件来克服. 第一种说法是, 就是我假想空间中有无数个相同的液体元胞的展开, 在立方体边界上的微粒, 依然可以受到邻近的元胞的作用. 第二种说法, 临近的元胞都是本体的镜像. 落实到实处就是, 一个微粒运动出了盒子, 将从盒子的另一边再穿进来. 有点像小时后玩的贪吃蛇吧.
    在这里插入图片描述

    如此一来, 盒子便没有了边界, 也没有了在表面的粒子, 这样的盒子很容易地用坐标的形式表现出来, 不用去储存真实的坐标. 但是我们仍有办法, 去得到这个元胞在真实体系中展开的情形. 折算到盒子内部坐标通常称之为wrap, 而按真实展开的叫做unwrap.
    有个非常关键的问题是, 这个非常小的无限循环的小东西, 能不能和宏观的系统有着相同的性质? 这取决于分子间作用力的范围和所考察的范围. 概括的说, 对于短程的强相互作用, 作用力距离不能超过盒子的一半. 这可能有一点抽象, 请看好Fig1.13, 如果粒子1的作用力距离超过一倍的盒子长, 那么会出现粒子自己对自己作用, 这显然是make sense的; 如果超过一半的盒长, 那么对于某些粒子, 粒子1在正方向将会对其作用一次, 从反方向同样会再次对它作用, make no science sense. 这个就是最小镜像原则(minimum image convention), 任意一个粒子的作用力仅被计算一次.

    对1来说, 5仅仅被计算了一次, 而不是两次对于长程作用力, 会有特殊的方法和技巧, 比如Ewald和PPPM, 去避开这些影响, 我们以后会谈到这个话题.

    不管前面说的怎么好, 经验告诉我们, 周期性边界条件还是会对热力学平衡和结构性质产生微弱的影响. 如果条件允许, 应该尽可能增加模拟的规模. 涉及到PBC的计算. 这里我们就不去考察rhombic dodecahedron这样复杂的元胞了, 就以立方体为例, 讲一下思路. 要注意的是, 计算坐标和计算距离是有所不同的.

    这里盒子是放在原点, 边长为L的立方体.计算坐标时, 我们不但想知道其在盒子中的位置, 也想知道实际上它跑到哪去了, 这个数据可以用来计算扩散等参数. 这时我们需要引入一个image flag的概念.

    这个值代表着原子穿出盒子的次数, 正方向穿出image flag+1, 反方向-1.当微粒从正方向跑出盒子, 那就意味着它跑进了下一个盒子, 坐标值将大于L. 此时应该减去L的值, 将其wrap进盒子, 同时image flag+1;

    如果从反方向跑了出去, 那就需要加一个L让他wrap进盒子, 同时image flag-1.计算距离时分为两种情况, 仅仅计算wrap体系中两个微粒间作用力, 那么就有, coord1-coord2.

    如果这个距离大于1/2L, 那L减去这个距离;

    如果小于-1/2L, 那就需要加上L. 总之很简单, 在纸上画个坐标轴取个特殊值就能搞明白.

    展开全文
  • 然而两个侧边的周期性边界条件不会处理,还望大神指点。 由此可见,此技术问题并不涉及固体或者是流体有限元问题,而是基础的周期性边界处理。 我最近采用的是这篇文献的处理方法 《周期性边界条件的一种处理方法...
  • tid=1304321&highlight=%E7%BC%96%E7%BB%87 ABAQUS 2018 推出了周期性边界官方插件直接调用 http://forum.simwe.com/forum.php?mod=viewthread&tid=1300570&...

    http://forum.simwe.com/forum.php?mod=viewthread&tid=1304321&highlight=%E7%BC%96%E7%BB%87

    ABAQUS 2018 推出了周期性边界官方插件  直接调用

     

    http://forum.simwe.com/forum.php?mod=viewthread&tid=1300570&highlight=%E7%BC%96%E7%BB%87----------------digimat 软件

     

    http://blog.sina.com.cn/s/blog_68d0921b0102v9yw.html ----abaqus动态分析的帖子

     

    请问有谁能讲讲如何定义复合材料的界面性能参数吗?最好是0厚度界面-----------http://forum.simwe.com/thread-1118126-1-1.html-------------- cohesive behavior可以

     

     

    请问如何设置really simple GUI  Dialog builder。想自己添加参数按键-------------------------------------ans:建立想要的对话框样式,在相应的keywords里输入程序里的关键字

     

    请教编织物纤维束间接触应如何施加-----------------http://forum.simwe.com/forum.php?mod=viewthread&tid=1283430&highlight=%E7%BC%96%E7%BB%87

    转载于:https://www.cnblogs.com/uncle7/p/10380935.html

    展开全文
  • 结果表明:在周期性边界条件下,围岩体各处的温度随时间周期性地波动,并且温度波动的振幅随围岩深度的延伸呈负指数规律变化,围岩体不同层面上的温度波具有振幅衰减和波动延迟的特点。综合以上分析,推导了巷道围岩体...
  • 程序边界条件

    2018-05-18 16:00:59
    ABAQUS提供了相当丰富的单元类型,材料属性等数据库可供用户选择,但是工程问题是千变万化的,为了满足用户的特殊工程要求,ABAQUS为用户提供了强大而又灵活的用户子程序接口(USER ...ABAQUS的一个子程序边界条件
  • 该算法是使用西奈半岛的台球模拟周期性边界条件下的Lorentz气体而构建的。 该程序模拟给定数量的点状粒子的轨迹,这些粒子可以自由地沿直线移动,直到与散射之一碰撞为止。 由于我们对粘度的传输过程感兴趣,因此...
  • 如果你喜欢 pdepe,但想解决周期性边界条件的问题,试试这个程序。 语法与pdepe略有不同,因此请查看两个示例文件以了解其用法。 详细信息:使用线和 ode15s 的方法(刚性求解器)。
  • 已更新完毕,本部分为原文第八章和第九章,初始设置问题以及边界条件定义两部分的翻译及讲解,边界条件定义部分,手册中写的很明白,但在不同的其他模型中使用并不完全相同,使用的时候不要教条。

    设置中的问题与定义边界条件

    Palabos用户文档 的第八章和第九章

    (一)Setting Up A Problem

    原始文档

    Attributing dynamics objects

    When constructing a new block-lattice, you must decide what type of collision is going to be executed on each of its cells, by assigning them a dynamics object. To avoid bugs related to cells without dynamics objects, the constructor of a block-lattice assigns a default-alue for the dynamics to all cells, the so-called background dynamics:

    // Construct a new block-lattice with a background-dynamics of type BGK.
    MultiBlockLattice3D<T,DESCRIPTOR> lattice( nx, ny, nz, new BGKdynamics<T,DESCRIPTOR>(omega) );
    

    After this, the dynamics of each cell can be redefined in order to adjust the behavior of the cells locally. The background dynamics differs from usual dynamics objects in an important way. To obtain memory savings, the constructor of the block-lattice creates only one instance of the dynamics object, and all cells refer to the same instance. If you modify a dynamics object on one cell, it changes everywhere. As an example, consider the relaxation parameter omega, which is stored in the dynamics, not in the cell. A function call like

    lattice.get(0,0,0).getDynamics().setOmega(newOmega); 
    

    would affect the value of the relaxation parameter on each cell of the lattice. Note that this particular behavior of the background dynamics, having several cells referring to the same object, cannot be reproduced by other means than constructing a new block-lattice. All other functions which override the background-dynamics on given cells with a new dynamics object are defined in such a way as to create an independent copy of the dynamics object for each cell.
    If the reference semantics of the background dynamics is contrary to your needs, you can re-define all dynamics objects right after constructing a block-lattice, to guarantee that all cells have an independent copy:

    // Override background-dynamics to guarantee and independent per-cell
    // copy of the dynamics object.
    defineDynamics( lattice, lattice.getBoundingBox(), new BGKdynamics<T,DESCRIPTOR> );
    

    There exist different variants of the function defineDynamics which you can use to adjust the dynamics of a group of cells (their syntax can be found in the Appendix Operations on the block-lattice):

    • One-cell version: Assign a new dynamics to just one cell.
      Example: tutorial/tutorial2/tutorial2_3.cpp.
    • BoxXD version: Assign a new dynamics to all cells within a rectangular domain.
      Example: showCases/multiComponent2d/rayleighTaylor2D.cpp, or the 3d example.
    • DotListXD version: Assign a new dynamics to several cells, listed individually in a dot-list structure.
      Example: codesByTopic/dotList/cylinder2d.cpp
    • Domain-functional version: Provide an analytical function which indicates the coordinates of cells which get a new dynamics object.
      Example: showCases/cylinder2d/cylinder2d.cpp
    • Bool-mask version: Specify the location of cells which get a new dynamics object through a Boolean mask, represented through a scalar-field.
      Example: codesByTopic/io/loadGeometry.cpp

    Initial values of density and velocity

    It is quite common to assign an initial value of velocity and density to a lattice by initializing all cells to an equilibrium distribution with the chosen density and velocity value. While this approach is not sufficient for time-dependent benchmark problems which depend critically on the initial state, it is most often sufficient to get a simulation started with reasonable initial values. The function initializeAtEquilibrium (see Appendix Operations on the blocklattice) is provided to initialize the cells within a rectangular-shaped domain at an equilibrium distribution. It comes in two flavors: one with a constant density and velocity within the domain (see the example examples/showCases/cavity2d or 3d), and one with a space-dependent value of these macroscopic values, specified through a userdefined function (see the example examples/showCases/poiseuille).
    To initialize cells in a different way, it is simplest to apply a custom operator to the cells with the help of one-cell functionals and indexed one-cell functionals introduced in Section Convenience wrappers for local operations.

    文档翻译

    设置动态对象

    当构造一个新的块格时,你必须通过分配一个动态对象来决定在每个单元格上执行什么类型的碰撞。为了避免与没有动态对象的单元格相关的bug,块格的构造器将动态的默认值赋给所有单元格,即所谓的后台动态:

    // 以BGK动力学类作为背景动力学构造一个新的块格。
    MultiBlockLattice3D<T,DESCRIPTOR> lattice( nx, ny, nz, new BGKdynamics<T,DESCRIPTOR>(omega) );
    

    在这之后,每个单元格的动态可以重新定义,以调整单元格的局部行为。背景动力学与一般的动力学对象有重要的区别。为了节省内存,块格的构造函数只创建一个动态类对象的实例,所有单元格都引用同一个实例。如果您在任一个单元格上修改动态类对象,它将到处更改。例如,考虑松弛参数,它存储在动态类对象中,而不是在单元格中。如以下函数调用

    lattice.get(0,0,0).getDynamics().setOmega(newOmega); 
    

    会影响晶格中每个单元的弛豫参数的值。请注意,这是背景动力学的特定行为,几个单元格指向同一个对象,不能通过其他方法重制,只能通过构造一个新的块晶格来实现。使用新的动态类对象覆盖给定单元上的背景动态的所有其他函数的定义方式是,为每个单元创建动态类对象的独立副本。
    如果背景动态的引用语义与你的需要相反,你可以在构建一个块格之后重新定义所有的动态对象,以保证所有的单元都有一个独立的副本:

    // 覆盖背景动态以保证每个单元都有动态对象的独立副本。
    defineDynamics( lattice, lattice.getBoundingBox(), new BGKdynamics<T,DESCRIPTOR> );
    

    你可以使用defineDynamics函数的不同变体来调整一组细胞的动态(它们的语法可以在block-lattice的附录操作中找到):

    • One-cell版本:将一个新的动态分配给一个细胞。
      例子:tutorial/tutorial2/tutorial2_3.cpp
    • BoxXD版本:将一个新的动态分配给矩形域中的所有单元格。
      例子:showCases/multiComponent2d/rayleighTaylor2D.cpp或 3d 案例
    • DotListXD版本:将一个新的动态分配给几个单元,以点列表结构单独列出。
      例子:codesByTopic/dotList/cylinder2d.cpp
    • Domain-functional版本:提供一个分析功能,指出细胞的坐标,得到一个新的动态对象。
      例子:showCases/cylinder2d/cylinder2d.cpp
    • Bool-mask版本:指定单元格的位置,这些单元格通过一个布尔掩码得到一个新的动态对象,通过一个标量场表示。
      例子:codesByTopic/io/loadGeometry.cpp
    密度和速度的初始值

    通过用选定的密度和速度值初始化所有的单元格使其达到平衡分布,为晶格分配速度和密度的初值是很常见的。虽然这种方法不适用于依赖于初始状态的时变基准测试问题,但它通常足以使仿真从合理的初始值开始。函数initializeAtEquilibrium(参见块格上的附录操作)被用来初始化矩形区域内的单元格,使其处于平衡分布。它有两种使用方式:一种是在域内具有恒定的密度和速度(参见示例examples/showCases/cavity2d或3d),另一种是通过用户定义的函数指定的这些宏观值的空间依赖性值(参见示例examples/showCases/poiseuille)。

    解释说明

    在构建块格时生成的动态类对象只有一个,在堆区分配对象空间,然后给块格中的每一个单元传递一个动态类对象的指针,所以只有一个对象,所有单元格共享。为特定单元格单独设置动态类,是通过在兑取新建一个相应的动态类对象并将其指针传递给对应单元格来覆盖只想背景动态类对象的指针,从而达到重新设置类对象的效果。One-cell、BoxXD、DotListXD、Domain-functional、Bool-mask都是指定更改所限定作用范围的参数,你可以把它们都想象成stl标准模板库中的iterator迭代器,不仅在动态类的设定中有使用,基本上所有其他的设置,都有他们的身影。
    密度和速度设置,整个块格的密度和速度的初始化,方腔流是将整个流区初始化为相同的密度和速度,也就是稳态模拟最有可能达到的平衡态附近,同理,泊肃叶流是将整个流区初始化为泊肃叶分布形式的密度和速度,也是对应的稳态模拟最可能达到的平衡态的附近,此种设置,都是为了能让模拟更快达到稳态结果,如果不这样设置,可能会导致模拟失败,非稳态的计算,我认为这个初始值并不是特别重要,主要的是在边界以及特定流动区域不要出现计算的死点,才是初始值和对应模拟参数需要关注的重点。

    (二)Defining Boundary Condition

    原始文档

    Overview

    The library Palabos currently implements a bunch of different boundary conditions for straight, grid-aligned boundaries. For general boundaries, the available options include stair-cased walls through bounce-back nodes, and higherorder accurate curved boundaries.
    IMPORTANT: Curved boundaries are available as of version 1.0, and are part of a fully-featured parallel stack, including parallel voxelization and I/O extensions. Curved boundaries are not yet documented. You can however find a full-featured example in showCases/aneurysm.

    Grid-aligned boundaries

    The class OnLatticeBoundaryConditionXD
    Some of the implemented boundary conditions are local (and are therefore implemented as dynamics classes), some are non-local (and are therefore implemented as data processors), and some have both local and nonlocal components. To offer a uniform interface to the various possibilities, Palabos offers the interface OnLatticeBoundaryConditionXD which is responsible for instantiating dynamics objects, adding data processors, or both, depending on the chosen boundary condition. The following line for example is used to chose a Zou/He boundary condition:

    OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition = createZouHeBoundaryCondition3D<T,DESCRIPTOR>();
    

    The four possibilities currently offered are (see palabos.org-bc for details):

    modelclass
    Regularized BCcreateLocalBoundaryConditionXD
    Skordos BCcreateInterpBoundaryConditionXD
    Zou/He BCcreateZouHeBoundaryConditionXD
    Inamuro BCcreateInamuroBoundaryConditionXD

    With each of these boundary conditions, it is possible to implement velocity Dirichlet boundaries (all velocity components have an imposed value) and pressure boundaries (the pressure is imposed, and the tangential velocity components are zero). Furthermore, various types of Neumann boundary conditions are proposed. In these case, a zero-gradient for a given variable is imposed with first-order accuracy by copying the value of this variable from a neighboring cell.

    Boundaries on a regular block domain
    The usage of these boundary conditions is somewhat tricky, because Palabos needs to know if a given wall is a straight wall or a corner (2D), or a flat wall, an edge or a corner (3D), and it needs to know the wall orientation. To simplify this, all OnLatticeBoundaryConditionXD objects have a method setVelocityConditionOnBlockBoundaries and a method setPressureConditionOnBlockBoundaries, to set up boundaries which are located on a rectangular domain. If all nodes on the boundary of a given 2D box boundaryBox implement a Dirichlet condition, use the following command:

    Box2D boundaryBox(x0,x1, y0,y1);
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, boundaryBox, locationOfBoundaryNodes );
    

    The first argument indicates the boxed shape on which the boundary is located, and the second argument specifies on which node a velocity condition is to be instantiated. If all boundary nodes are located on the exterior bounding box of the lattice, the above command simply becomes:

    boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice, locationOfBoundaryNodes);
    

    Finally, if simply all nodes of the exterior bounding box shall implement a velocity condition, this simplifies to:

    boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice);
    

    Now, suppose that the upper and lower walls, including corners, in a nx -by- ny lattice implement a free-slip condition (zero-gradient for tangential velocity components), the left wall a Dirichlet condition, and the right wall and outflow condition (zero-gradient for all velocity components). This is achieved by adding one more parameter to the function call:

    Box2D inlet(0, 0, 2, ny-2);
    Box2D outlet(nx-1, nx-1, 2, ny-2);
    Box2D bottomWall(0, nx-1, 0, 0);
    Box2D topWall(0, nx-1, ny-1, ny-1);
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, inlet );
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, outlet, boundary::outflow );
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, bottomWall, boundary::freeslip );
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, topWall, boundary::freeslip );
    

    Remember that if the boundary is not located on the outer bounds of the lattice, you must use an additional Box3D argument to say where the boundaries are located, additionally to saying which boundary is going to be added. This is necessary, because Palabos needs to know the orientation and the nature (flat wall, corner, etc.) of the boundary nodes.

    Zero-gradient conditions
    The optional arguments like boundary::outflow are of type boundary::BcType, and can have the following values for velocity boundaries:

    object valuesexplain
    boundary::dirichlet (default value)Dirichlet condition: imposed velocity value
    boundary::outflow or boundary::neumannZero-gradient for all velocity components
    boundary::freeslipZero-gradient for tangential velocity components, zero value for the others
    boundary::normalOutflowZero-gradient for normal velocity components, zero value for the others

    For pressure boundaries, you have the following choice:

    object valuesexplain
    boundary::dirichlet (default value)Dirichlet condition: imposed pressure value, zero value for the tangential velocity components
    boundary::neumannZero-gradient for the pressure, zero value for the tangential velocity components

    Boundary conditions cannot be overriden
    It is not possible to override the type of a boundary. Once a boundary node is a Dirichlet velocity node, it stays a Dirichlet velocity node. The result of re-defining it as, say, a pressure boundary node, is undefined. Therefore, boundary conditions must be defined carefully, and piece-wise, as shown in the example above. On the other hand, the value imposed on the boundary (i.e. the velocity value on a Dirichlet velocity boundary) can be changed as often as needed.

    Setting the velocity or pressure value on Dirichlet boundaries
    A constant velocity value is imposed with the function setBoundaryVelocity. The following line sets the velocity values to zero on all nodes of a 2D domain which have previously been defined to be Dirichlet velocity nodes:

    setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array<T,2>(0.,0.) );
    

    On all other nodes, this command has no effect. To set a constant pressure value (or equivalently, density value), use the command setBoundaryDensity:

    setBoundaryDensity(lattice, lattice.getBoundingBox(), 1. );
    

    A non-constant, space dependent velocity resp. pressure profile can be defined by replacing the velocity resp. pressure argument by either a function or by a function object (an instance of a class which overrides the function call operator) which yields a value of the velocity resp. density as a function of space position. An example is provided in the file examples/showCases/poiseuille/poiseuille.cpp.

    Interior and exterior boundaries
    Sometimes, the boundary is not located on the surface of a regular, block-shaped domain. In this case, the functions like setVelocityConditionOnBlockBoundaries are useless, and the boundary must be manually constructed. As an example, here’s how to construct manually a free-slip condition along the top-wall, including the corner nodes:

    // The automatic creation of a top-wall ...
    Box2D topWall(0, nx-1, ny-1, ny-1);
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, topWall, boundary::freeslip );
    
    // ... can be replaced by a manual construction as follows:
    Box2D topLid(1, nx-2, ny-1, ny-1);
    boundaryCondition->addVelocityBoundary1P( topLid, lattice, boundary::freeslip );
    boundaryCondition->addExternalVelocityCornerNP( 0, ny-1, lattice, boundary::freeslip );
    boundaryCondition->addExternalVelocityCornerPP( nx-1, ny-1, lattice, boundary::freeslip );
    

    A distinction is made between external and internal corners, depending on whether the corner is convex or concave. On the following example geometry, you’ll find five external and one internal corner:
    the picture of boundary direction
    The extensions like 1P, NP, and PP at the end of the methods of the boundary-condition object are used to indicate the orientation of the wall normal, pointing outside the fluid domain, as shown on the figure above. On a straight wall, the code 1P means: “the wall normal points into positive y-direction”. Likewise, the inlet would be labeled with the code 0N as in “negative x-direction”. On a corner, the code NP means “negative x-direction and positive y-direction”. It is important to mention that this wall normal is purely geometrical and does not depend on whether a given wall has the function of an inlet or an outlet. In both cases, the wall normal points away from the fluid, into the non-fluid area.
    In 3D, you’ll use the following function for plane walls:

    addVelocityBoundaryDO
    

    where the direction D can be 0 for x, 1 for y, or 2 for z, and the orientation O has the value P for positive and N for negative. Edges can be internal or external. For example:

    addInternalVelocityEdge0NP
    addExternalVelocityEdge1PN
    

    where the first digit of the code indicates the axis to which the edge is parallel, and the two subsequent digits indicate the orientation of the edge inside the plane normal to the edge, in the same way as the corner nodes in 2D. The axes are counted periodically: 0NP means x-plane, negative y-value, and positive z-value, whereas 1PN means y-plane, positive z-value, and negative x-value. Finally, 3D corners are constructed in the same way as in 2D, through a function call like the following:

    addInternalVelocityCornerPNP
    addExternalVelocityCornerNNN
    

    Periodic boundaries

    The concept of periodicity in Palabos is different from the approach chosen for the other types of boundary conditions. Periodicity can only be implemented on opposite, outer boundaries of an atomic-block or a multi-block. This behavior works also if the multi-block has a sparse memory implementation. In this case, all fluid nodes which are in contact with an outer boundary of the multi-block can be periodic, if the corresponding node on the opposite wall is a fluid node.
    In the case of an atomic-block (remember that this case is uninteresting, because you should always work with multiblocks anyway), periodicity is simply a property of the streaming operator. It has no effect for scalar-fields and tensor-fields. In the case of a multi-block however, periodicity can also affect scalar-fields and tensor-fields. Being periodic in this case means that if you define a non-local data processor which accesses nearest neighbor nodes, the value of these neighbor nodes is determined by periodicity along an outer boundary of the multi-block.
    By default, all blocks are non-periodic: as mentioned previously, the behavior of outer boundary nodes defaults to one of the versions of bounce-back algorithm. To turn on periodicity along, say, the x-direction, write a command like

    lattice.periodicity().toggle(0, true); // Periodic block-lattice.
    scalarField.periodicity().toggle(0, true); // Periodic scalar-field.
    

    You can also define all boundaries to be periodic through a single command:

    lattice.periodicity().toggleAll(true);
    

    Please note that the periodicity or non-periodicity of a block should be defined as soon as possible after the creation of the block, because Palabos needs to know if boundaries are periodic or not when adding or executing data processors.
    As a last remark, you should be aware that boundaries of Neumann type (outflow, free-slip, adiabatic thermal wall, etc.) should never be defined to be periodic. This wouldn’t make sense anyway, and it leads to an undefined behavior (aka program crash).

    Bounce-back

    In Palabos, all outer boundaries of a lattice which are not periodic automatically implement a version of the bounceback boundary algorithm, according to the following algorithm. In pre-collision state, all unknown populations are assigned the post-collision value on the same node at the end of the previous iteration, from the opposite location. We will call this algorithm Bounce-Back 1.
    A bounce-back condition can also be used elsewhere in the domain, by explicitly assigning a dynamics object of type BounceBack<T,Descriptor> to chosen cells. On these newly assigned nodes, the collision step is replaced by a bounce-back procedure, during which each population is replaced by the population corresponding to the opposite direction. Such a node can not be considered any more as a fluid node. Computing velocity-moments of the populations on a bounce-back node yields for example arbitrary results, because some of the populations (those which are not incoming from the fluid) have arbitrary values. Consequently, you cannot compute macroscopic quantities like density and velocity on these nodes. Instead, these nodes should be considered as pure algorithmic entities which have the purpose to re-inject into the fluid all populations that leave the domain. Consider fluid nodes adjacent to this type of bounce-back cells. If you think about it, you’ll realize that these nodes behave exactly as if they were cells of type Bounce-Back 1, with a small difference: unknown populations at time t get the post-collision value from opposite populations from time t-2 (and not time t-1 as in Bounce-Back 1. We will refer to this version of bounce-back as Bounce-Back 2.
    The effect of both versions of bounce-back is to produce a no-slip wall located half-way between nodes. In the case of Bounce-Back 1, this is half a cell-width beyond the bounce-back fluid cell, whereas in the case of Bounce-Back 2 it is half-way between the last fluid cell and the non-fluid bounce-back node. Obviously, Bounce-Back 2 wastes numerical precision over Bounce-Back 1 by leaping over a time step, and therefore leads to a lower-order accuracy in time of the numerical method. Furthermore, both versions of bounce-back represent the wall location through a staircase approximation, which is low-order accurate (“Order-h”) in space. On the other hand, Bounce-Back 2 offers a huge advantage: its implementation is independent of the wall orientation. If you decide that a node is a bounce-back node, you simply say so; no need to know if it is a corner, a left wall, or a right wall. Constructing a geometry is then as easy as filling areas with bounce-back nodes. In many cases, this ease of use makes up, to a large extent, for the low-order accuracy of the wall representation. When setting up a new simulation, it is always good to try Bounce-Back 2 as a first attempt, as the quality of the obtained results is often surprising. In highly complex geometry such as porous media, bounce-back is often even considered superior to other approaches.

    Bounce-back domain from analytical description
    The Rule 0 in Palabos can be stated often enough: don’t write loops over space in end-user code. This also true, of course, when assigning bounce-back dynamics to lattice cells. To guarantee reasonable performance, this is task is performed inside a data processor, or even better, by using the wrapper function defineDynamics. This process is illustrated in the example program examples/codesByTopic/bounceBack/instantiateCylinder.cpp. First, a class is written which defines the location of the bounce-back nodes as a function of the cell coordinates iX and iY. In the present case, the nodes are located inside a 2D circular domain:

    template<typename T>
    class CylinderShapeDomain2D : public DomainFunctional2D {
    public:
    	CylinderShapeDomain2D(plint cx_, plint cy_, plint radius) 
    						: cx(cx_), cy(cy_), radiusSqr(util::sqr(radius)) { }
    	// The function-call operator is overridden to specify the location
    	// of bounce-back nodes.
    	virtual bool operator() (plint iX, plint iY) const {
    		return util::sqr(iX-cx) + util::sqr(iY-cy) <= radiusSqr;
    	}
    	virtual CylinderShapeDomain2D<T>* clone() const {
    		return new CylinderShapeDomain2D<T>(*this);
    	}
    private:
    	plint cx, cy;
    	plint radiusSqr;
    };
    

    An instance of this class is then provided as an argument to defineDynamics:

    defineDynamics( lattice, lattice.getBoundingBox(),
    				new CylinderShapeDomain2D<T>(cx, cy, radius),
    				new BounceBack<T,DESCRIPTOR> );
    

    The second argument is of type Box2D, and it can be used to improve the efficiency of this function call: the bounce-back nodes are attributed on cells on the intersection between this box and the domain specified by the domain-functional.
    If your domain is specified as the union of a collection of simpler domains, you can use defineDynamics iteratively on each of the simpler domains. If the domain is the intersection or any other geometric operation of simpler domains, you’ll need to play with boolean operators inside the definition of your domain-functional.

    Bounce-back domain from boolean mask
    In this section, it is assumed that the geometry of the domain is prescribed from an external source. This is for example the case when you simulate a porous media and possess data from a scan of the porous material in question. This easiest approach consists in storing this data as an ASCII file which distinguishes fluid nodes from solid nodes by zeros and ones, separated with spaces, tabs, or newline characters. The file represents the content of a regular array and stores only raw data, no information on the matrix dimensions. The data layout must be the same as in Palabos: an increment of the z-coordinate represents a continuous progress in the file (in 3D), followed by the y-coordinate, and then the x-coordinate. This data is simply flushed from the file into a scalar-field. The following code is taken from the example examples/codesByTopic/io/loadGeometry.cpp:

    MultiScalarField2D<T> boolMask(nx, ny);
    plb_ifstream ifile("geometry.dat");
    ifile >> boolMask;
    

    Note that currently, no error-checking is implemented for such I/O operations. It is your responsibility to ensure that the dimensions of the scalar-field correspond to the size of the data in the geometry file. As a next step, the bounce-back nodes are instantiated using one of the versions of the function defineDynamics:

    defineDynamics(lattice, boolMask, new BounceBack<T,DESCRIPTOR>, true);
    

    This function is currently only defined for multi-blocks using the same data type (T in this case), which explains why the scalar-field boolMask is of type T and not bool, as one could have expected. The scalar-field is not bool, as one could have expected, but the real type T.

    文档翻译

    综述

    Palabos库目前为网格对齐的直线边界实现了许多不同的边界条件。对于一般的边界,可用的选项包括通过反弹节点的阶梯式墙壁,以及更高阶的精确弯曲边界。
    重要提示:从1.0版开始就可以使用曲线边界,它是一个功能齐全的并行堆栈的一部分,包括并行体素化和I/O扩展。弯曲的边界还没有记录。然而,你可以在showCases/aneurysm中找到一个全功能的例子。

    网格对齐边界

    OnLatticeBoundaryConditionXD类
    一些实现的边界条件是本地的(因此被实现为动态类),一些是非本地的(因此被实现为数据处理器),一些同时具有本地和非本地组件。为了提供各种可能性的统一接口,Palabos提供OnLatticeBoundaryConditionXD接口,该接口负责实例化动态对象、添加数据处理器,或两者兼得,这取决于所选择的边界条件。以下面这行为例,选择一个邹/河边界条件:

    OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition = createZouHeBoundaryCondition3D<T,DESCRIPTOR>();
    

    目前提供的四种可能性是:

    modelclass
    Regularized BCcreateLocalBoundaryConditionXD
    Skordos BCcreateInterpBoundaryConditionXD
    Zou/He BCcreateZouHeBoundaryConditionXD
    Inamuro BCcreateInamuroBoundaryConditionXD

    对于每一个边界条件,都可以实现速度Dirichlet边界(所有速度分量都有一个施加的值)和压力边界(施加压力,切向速度分量为零)。此外,还提出了各种类型的Neumann边界条件。在这种情况下,通过从邻近单元格复制该变量的值,以一阶精度对给定变量施加零梯度。

    常规块域上的边界
    这些边界条件的使用有些棘手,因为Palabos需要知道给定的墙是直墙还是角(2D),还是平墙、边或角(3D),并且需要知道墙的方向。为了简化这一点,所有OnLatticeBoundaryConditionXD对象都有一个方法setVelocityConditionOnBlockBoundaries和一个方法setPressureConditionOnBlockBoundaries,用来设置位于矩形域中的边界。如果给定二维框边界上的所有节点都实现Dirichlet条件,则使用以下命令:

    Box2D boundaryBox(x0,x1, y0,y1);
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, boundaryBox, locationOfBoundaryNodes );
    

    第一个参数表示边界所在的方箱形状,第二个参数指定要在哪个节点上实例化velocity条件。如果所有边界节点都位于格外边界框内,则上述命令简单为:

    boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice, locationOfBoundaryNodes);
    

    最后,如果外边界框的所有节点都要实现速度条件,则简化为:

    boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice);
    

    现在,假设在一个二维晶格中,上下壁面(包括角)实现了一个自由滑动条件(切向速度分量为零梯度),左壁面为一个Dirichlet条件,右壁面和出流条件(所有速度分量为零梯度)。这是通过在函数调用中增加一个参数来实现的:

    Box2D inlet(0, 0, 2, ny-2);
    Box2D outlet(nx-1, nx-1, 2, ny-2);
    Box2D bottomWall(0, nx-1, 0, 0);
    Box2D topWall(0, nx-1, ny-1, ny-1);
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, inlet );
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, outlet, boundary::outflow );
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, bottomWall, boundary::freeslip );
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, topWall, boundary::freeslip );
    

    请记住,如果边界不在格的外部边界上,则必须使用额外的Box3D参数来指定边界的位置,此外还要说明将添加哪个边界。这是必要的,因为Palabos需要知道边界节点的方向和性质(平墙、角等)。

    零梯度条件
    The optional arguments like boundary::outflow are of type boundary::BcType, and can have the following values for velocity boundaries:
    可选参数如boundary::outflow的类型为boundary::BcType,速度边界可以有以下值:

    object valuesexplain
    boundary::dirichlet (default value)Dirichlet 条件: 实施速度值
    boundary::outflow or boundary::neumann所有速度分量零梯度
    boundary::freeslip切向速度分量为零梯度,其它分量为零
    boundary::normalOutflow法向速度分量为零梯度,其他分量为零

    对于压力边界,你有以下选择:

    object valuesexplain
    boundary::dirichlet (default value)Dirichlet 条件: 施加压力值,切向速度分量为零
    boundary::neumann压强为零梯度,切向速度分量为零

    边界条件不能被覆盖
    不可能重写边界的类型。一旦边界节点是狄利克雷速度节点,它仍然是狄利克雷速度节点。将其重新定义为,例如,一个压力边界节点的结果是未定义的。因此,边界条件必须谨慎地、分段地定义,如上例所示。另一方面,施加在边界上的值(即狄利克雷速度边界上的速度值)可以根据需要随时改变。

    在狄利克雷边界上设置速度或压力值
    由函数setBoundaryVelocity施加一个恒定的速度值。下面的线将二维域中所有节点上的速度值设为0,之前定义为狄利克雷速度节点:

    setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array<T,2>(0.,0.) );
    

    在所有其他节点上,此命令无效。要设置一个恒定的压力值(或等效的密度值),使用命令setBoundaryDensity:

    setBoundaryDensity(lattice, lattice.getBoundingBox(), 1. );
    

    A non-constant, space dependent velocity pressure profile can be defined by replacing the velocity resp. pressure argument by either a function or by a function object (an instance of a class which overrides the function call operator) which yields a value of the velocity resp. density as a function of space position. An example is provided in the file examples/showCases/poiseuille/poiseuille.cpp.
    一个非恒定的,依赖于空间的速度分布和压力分布可以分别通过替换速度分布和压力分布来确定。一个函数或一个函数对象(一个覆盖函数调用操作符的类的实例)的压力参数,该参数分别产生速度和的值作为空间位置的函数。示例在examples/showCases/poiseuille/poiseuile .cpp文件中提供。

    内外边界
    有时,边界并不位于规则的块状区域的表面。在这种情况下,像setVelocityConditionOnBlockBoundaries这样的函数是没有用的,必须手工构造边界。作为一个例子,这里展示如何手动构建一个沿顶壁自由滑动的边界条件,包括角节点:

    // The automatic creation of a top-wall ...
    Box2D topWall(0, nx-1, ny-1, ny-1);
    boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, topWall, boundary::freeslip );
    
    // ... can be replaced by a manual construction as follows:
    Box2D topLid(1, nx-2, ny-1, ny-1);
    boundaryCondition->addVelocityBoundary1P( topLid, lattice, boundary::freeslip );
    boundaryCondition->addExternalVelocityCornerNP( 0, ny-1, lattice, boundary::freeslip );
    boundaryCondition->addExternalVelocityCornerPP( nx-1, ny-1, lattice, boundary::freeslip );
    

    外部和内部的角是有区别的,这取决于角是凸的还是凹的。在下面的几何示例中,你会发现五个外部和一个内部角落:
    the picture of boundary direction
    如上图所示,边界条件对象方法末尾的1P、NP和PP等扩展用于指示壁面法线指向流体域外的方向。在一面直墙上,代码1P表示:“墙法线指向正的y方向”。同样地,入口将被标记代码0N为“负x方向”。在一个角落里,NP代表“负的x方向和正的y方向”。需要指出的是,这道墙法线完全是几何图形,并不取决于给定的墙是具有入口还是出口的功能。在这两种情况下,壁法线点远离流体,进入非流体区。
    在3D中,你将使用以下功能来制作平面墙:

    addVelocityBoundaryDO
    

    方向D可以是x为0,y为1,z为2,方向O为P为正,N为负。边可以是内部的,也可以是外部的。例如:

    addInternalVelocityEdge0NP
    addExternalVelocityEdge1PN
    

    其中,代码的第一个数字表示与边缘平行的轴,随后的两个数字表示与边缘垂直的平面内边缘的方向,与2D中的角节点相同。轴是周期性计数的:0NP表示x平面,负的y值,正的z值,而1PN表示y平面,正的z值,负的x值。最后,通过如下所示的函数调用,以与2D相同的方式构造3D角:

    addInternalVelocityCornerPNP
    addExternalVelocityCornerNNN
    
    周期边界

    Palabos中周期性的概念不同于其他类型边界条件所选择的方法。周期性只能在原子块或多块的相对、外部边界上实现。如果多块具有稀疏内存实现,则此行为也有效。在这种情况下,如果边界对面壁上的对应节点是流体节点,则与多块外边界接触的所有流体节点都是周期性的。
    在原子块的情况下(请记住,这种情况是无趣的,因为无论如何都应该使用多块),周期性只是流操作符的一个属性。它对标量场和张量场没有影响。然而,在多块的情况下,周期性也会影响标量域和张量域。在这种情况下,周期性意味着如果定义一个访问最近邻居节点的非本地数据处理器,这些邻居节点的值将由沿多块外部边界的周期性决定。
    默认情况下,所有块都是非周期性的:如前所述,外部边界节点的行为默认为一个反弹算法的版本。要打开沿x方向的周期性,可以这样编写命令

    lattice.periodicity().toggle(0, true); // Periodic block-lattice.
    scalarField.periodicity().toggle(0, true); // Periodic scalar-field.
    

    你也可以定义所有的边界是周期性的通过一个单一的命令:

    lattice.periodicity().toggleAll(true);
    

    请注意,块的周期性或非周期性应该在块创建后尽快定义,因为Palabos需要知道添加或执行数据处理器时边界是否是周期性的。
    最后,你应该知道Neumann类型的边界(流出、自由滑动、绝热热壁等)不应该被定义为周期性的。无论如何,这都没有意义,并且会导致未定义的行为(即程序崩溃)。

    反弹

    在Palabos中,晶格中所有非周期的外边界都会根据下面的算法自动实现一个版本的反弹边界算法。在碰撞前状态下,所有未知种群在前一次迭代结束时,从相反的位置在同一节点上分配碰撞后值。我们称这个算法为Bounce-Back 1。
    通过显式地将类型为BounceBack<T,Descriptor> 分配给所选的单元格,还可以在域中的其他地方使用反弹条件。在这些新分配的节点上,碰撞步骤将被一个反弹过程替换,在此过程中,每个种群将被对应于相反方向的种群替换。这样的节点不能再被认为是流动节点。例如,在回弹节点上计算种群的速度矩可以得到任意的结果,因为有些种群(那些不是来自流体的种群)具有任意的值。因此,你不能在这些节点上计算像密度和速度这样的宏观量。相反,这些节点应该被视为纯粹的算法实体,其目的是将所有离开域的种群重新注入流体。考虑与这种弹回单元相邻的流体节点。如果你仔细想想,你会发现这些节点的行为与Bounce-Back 1类型的单元完全一样,只是有一点不同:t时刻的未知种群从t-2时刻的相反种群获得碰撞后的值(而不是像弹回1那样从t-1时刻获得碰撞后的值)。我们将把这个版本的反弹称为Bounce-Back 2。
    这两个版本的弹回的效果是产生一个位于中间的无滑动的墙节点。在Bounce-Back 1的情况下,这是反弹流体单元的一半宽度,而在Bounce-Back 2的情况下,这是在最后一个流体单元和非流体反弹节点之间的一半宽度。很明显,由于跨越了时间步长,Bounce-Back 2比Bounce-Back 1浪费了数值精度,因此导致数值方法在时间上的低阶精度。此外,这两个版本的反弹代表了通过楼梯近似值的墙壁位置,这是低阶精度(“Order-h”)在空间。另一方面,Bounce-Back 2提供了一个巨大的优势:它的实现独立于墙的方向。如果您确定某个节点是回跳节点,可以简单的这样说;不需要知道它是一个角落,左边的墙,还是右边的墙。然后,构造一个几何图形就像用反弹节点填充区域一样简单。在许多情况下,这种易用性在很大程度上弥补了墙壁表示的低阶精度。在设置新的模拟时,最好先尝试一下Bounce-Back 2,因为得到的结果的质量常常令人惊讶。在高度复杂的几何结构中,如多孔介质,反弹甚至常常被认为优于其他方法。

    解析描述反弹域
    Palabos中的规则0需要经常声明:不要在最终用户代码中编写空间上的循环。当然,当将反弹动力学分配给晶格单元时,也是如此。为了保证合理的性能,这是在数据处理器中执行的任务,或者更好,使用包装器函数defineDynamics。这个过程在示例程序examples/codesByTopic/bounceBack/instantiateCylinder.cpp中进行了说明。首先,编写一个类,该类将反弹节点的位置定义为单元坐标iX和iY的函数。在本例中,节点位于二维圆形域内:

    template<typename T>
    class CylinderShapeDomain2D : public DomainFunctional2D {
    public:
    	CylinderShapeDomain2D(plint cx_, plint cy_, plint radius) 
    						: cx(cx_), cy(cy_), radiusSqr(util::sqr(radius)) { }
    	// 函数调用操作重写被指定反弹节点的位置。
    	virtual bool operator() (plint iX, plint iY) const {
    		return util::sqr(iX-cx) + util::sqr(iY-cy) <= radiusSqr;
    	}
    	virtual CylinderShapeDomain2D<T>* clone() const {
    		return new CylinderShapeDomain2D<T>(*this);
    	}
    private:
    	plint cx, cy;
    	plint radiusSqr;
    };
    

    然后,这个类的一个实例作为参数提供给defineDynamics:

    defineDynamics( lattice, lattice.getBoundingBox(),
    				new CylinderShapeDomain2D<T>(cx, cy, radius),
    				new BounceBack<T,DESCRIPTOR> );
    

    第二个参数是Box2D类型的,它可以用来提高这个函数调用的效率:在这个框和域函数指定的域之间的交集上的单元格上赋予反弹节点属性。
    如果您的域被指定为一组更简单的域的并集,那么您可以在每个更简单的域上迭代地使用defineDynamics。如果域是简单域的交集或其他几何运算,则需要在域函数的定义中使用布尔运算符。

    从布尔掩码返回域
    在本节中,假设定义域的几何形状是由外部源规定的。例如,当您模拟多孔介质并拥有所述多孔材料的扫描数据时,就会出现这种情况。这种最简单的方法是将数据存储为ASCII文件,该文件通过0和1区分流动节点和实体节点,用空格、制表符或换行符分隔。该文件表示常规数组的内容,只存储原始数据,不存储关于矩阵维的信息。数据布局必须与Palabos相同:z坐标的增量表示文件中的连续进度(在3D中),接着是y坐标,然后是x坐标。这些数据只是从文件中刷新到一个标量字段中。下面的代码摘自例子 examples/codesByTopic/io/loadGeometry.cpp:

    MultiScalarField2D<T> boolMask(nx, ny);
    plb_ifstream ifile("geometry.dat");
    ifile >> boolMask;
    

    请注意,目前没有为此类I/O操作实现错误检查。您有责任确保标量字段的大小与几何文件中的数据大小一致。下一步,使用函数defineDynamics的一个版本实例化回跳节点:

    defineDynamics(lattice, boolMask, new BounceBack<T,DESCRIPTOR>, true);
    

    这个函数目前仅使用相同的数据类型(在本例中为T)为多个块定义,这解释了为什么scalar-field布尔掩码的类型是T而不是bool,正如人们所期望的那样。标量场不是人们所期望的bool,而是真正的T类型。

    解释说明

    对于边界条件,操作方法都在这了,基本没什么遗漏,反弹边界条件是通过动态类设置的,是可以覆盖设置的,OnLatticeBoundaryConditionXD类,如果没看错,是通过使用枚举值调用特定的函数来实现的,不能重复设置,请注意,一般情况下,因为我主要是做多孔介质方面,一个一个边角点设置边界没有可能,我选用的就是最基本的反弹边界条件。
    因为并没有做到,我并没有继续翻译后面一小段的曲面边界的内容。如果有需要,请自行观看实验。
    还应注意,不同的模型对于边界的处理也是不同的,相应的调用的枚举值也是不一样的,需要各位在仔细阅读相关的模型案例后,自行判断使用方法。

    展开全文
  • 对商业周期的分析需要提取时间序列的周期性成分,该时间序列通常也受到诸如潜在趋势或噪声等其他因素的影响。本文介绍了一些在最近的文献中用于从给定系列中提取商业周期的方法。它基于Stock and Watson(1999)在...

    原文链接:http://tecdat.cn/?p=5399

    原文出处:拓端数据部落公众号

    介绍

    对商业周期的分析需要提取时间序列的周期性成分,该时间序列通常也受到诸如潜在趋势或噪声等其他因素的影响。本文介绍了一些在最近的文献中用于从给定序列中提取商业周期的方法。它基于Stock and Watson(1999)在“宏观经济学手册”中关于商业周期的章节。我还介绍了相对较新的方法,如小波滤波器或经验模式分解。由于这篇文章的重点是在R中实现某些过滤技术,我不会涉及数学。相反,我将参考各自的文献。对于这些例子,我使用了实际GDP的季度数据。

    names(gdp) <- c("Time","GDP") # 重命名
    
    gdp[,"GDP"] <- log(gdp[,"GDP"]) # 取对数

    为了直观地了解提取时间序列的周期性成分意味着什么,请查看下图中随时间变化的对数实际GDP的发展情况。

    
    ggplot(gdp,aes(x=Time,y=GDP)) + geom_line(size=.5) + theme_classic() 

    数据有明显的增长趋势,逐渐变小。此外,该序列以一种或多或少的常规方式围绕这一趋势波动。该序列与趋势的偏差非常小,这种偏差经常发生,但也有相当大的偏差,这种偏差可能会持续几个时期。后者是与商业周期分析相关的波动。

    时间趋于衰退

    从一序列中剔除趋势的第一种方法是在时间变量上回归感兴趣的变量并获得残差。这些在下图中绘制,其中线性趋势被移除。

    # 可视化
    
    
    ggplot(dat,aes(x=Time,y=Linearly.Detrended)) + geom_hline+ geom_line

    这种方法相对有争议,因为它假设存在一个恒定的线性时间趋势。正如我们上面所看到的,鉴于趋势增长率随时间的稳步下降,这种情况不太可能发生。然而,仍然可以假设时间趋势的不同函数形式,例如添加二次项,摆脱趋势。这种方法的另一个缺点是它只能排除趋势,而不能排除噪声,即序列中的非常小的波动。

    差分

    接下来的方法是采用一阶差分,因为它通常被教导以获得平稳的时间序列。这假设数据是不稳定的。取得一阶差分的结果显示在下图中,其中它也与时间趋势序列进行比较。差分数据在零线附近波动得更多,但它也包含很多噪声。

    #可视化
    
     g <- melt(dat,,na.rm=TRUE)
    
    
    
    
    # 定义绘图函数
    
    plot
    
    
    axis.line=element_line(size=.3,colour="black"), # 设置绘图参数
    
    
    
    
    # 绘图
    
    plot.cycles(d=g)

    Hodrick Prescott过滤器

    Hodrick和Prescott(1981)开发了一种滤波器,它将时间序列分为趋势,周期和噪声分量。需要时间序列和平滑参数。下图显示了Hodrick-Prescott滤波器获得的实际GDP的周期性成分值,并将其与线性去趋势序列的值进行了比较。两个序列的行为看起来非常相似,只是HP序列在零附近波动较大,而线性去趋势序列仍然包含趋势的成分。此外,循环HP序列还包括一些类似噪音的成分

    #绘图
    
    
    g <- melt(dat[,c(1,4,3)],id.vars="Time",na.rm=TRUE)
    
    
    plot(g)

    Baxter过滤器

    Baxter和King(1994,1999)提出了一种滤波器,它可以生成与HP滤波器类似的结果,但它可以剔除上面显示的许多类似噪声的行为。它需要序列,周期数量的下限和上限,假定周期发生(plpu),以及平滑因子nfix。文献(参见NBER,Stock和Watson(1999))表明商业周期持续6至32个月。这些值用于指定循环周期的下限和上限。BK滤波器的结果如下图所示。该方法的一个相对序列的缺点是平滑因子导致在序列的开始和结束时样本的丢失。这可能是小样本的问题。

    # 绘图
    
    dat <- cbind(dat,data.frame( bk))
    
    g <- melt(dat[,c(1,5,4)], ,na.rm=TRUE)
     
    
    plot(g )

    小波滤波器

    Yogo(2008)提出使用小波滤波器从时间序列数据中提取商业周期。这种方法的优点是该功能不仅可以提取序列的趋势,周期和噪声,而且可以更加具体地说明周期发生的周期。然而,由于该方法只能捕获2的幂的周期性,即2,4,8,16,32等,所以没有完全的自由度。

    R中的方法实现也很简洁,但在使用之前需要一些额外的数据转换。它需要时间序列的不同版本和分解的深度。

    该函数给出了多个序列,必须将它们累积起来cumsum,将它们转换回反映周期性模式的数据。此外,一些序列可以结合使用rowSums。当应该一起分析持续8到16和16到32个周期的周期时,这很有用,如下图所示。小波滤波器产生与BK滤波器类似的结果,因为循环周期的上限在两者中相等,下限仅相差2。

    #绘图
    
     g <- melt(dat[,c(1,6,5)],id.vars="Time",na.rm=TRUE)
    
    levels(g[,2]) <- c("Wavelet","Baxter King")
    
    
    
    plot.cycles(g,"Wavelet vs. Baxter King")

    经验模式分解(EMD)

    基于Huang等人。(1998)Kozic和Sever(2014)提出经验模式分解作为商业周期提取的另一种方法。它需要不同的时间序列,边界条件和规则,该规则指定迭代过程在哪个点获得了满意的结果并且可以停止。该滤波器方法的结果与HP,BK和小波滤波器相比有所不同。

     emd <- as.data.frame(emd(xt=diff(gdp[,2]),boundary="wave",stoprule="type2")$imf)
    
    
    
     g <- melt(dat[,c(1,7,4)],id.vars="Time",na.rm=TRUE)
    
      
    
    plot.cycles(g,"EMD vs. Hodrick Prescott")

     


    最受欢迎的见解

    1.在python中使用lstm和pytorch进行时间序列预测

    2.python中利用长短期记忆模型lstm进行时间序列预测分析

    3.使用r语言进行时间序列(arima,指数平滑)分析

    4.r语言多元copula-garch-模型时间序列预测

    5.r语言copulas和金融时间序列案例

    6.使用r语言随机波动模型sv处理时间序列中的随机波动

    7.r语言时间序列tar阈值自回归模型

    8.r语言k-shape时间序列聚类方法对股票价格时间序列聚类

    9.python3用arima模型进行时间序列预测

    展开全文
  • 没多久我就决定要开始我的刷题之路,我仍然记得第一题好像是求两数之和,其实这个题目我在备考时的九度OJ上做过,感觉应该难度不算太大,磨了一些边界之后,当一个绿色的Accept出现在眼前,我觉得这种感觉就像是我...
  • 为形成一套完善的轴棒法炭/炭复合材料基本力学性能预测方法,首先完成了轴棒法炭/炭复合材料代表性体积单元的建立,然后以周期性边界条件为理论基础,编写了周期性边界条件程序,最后结合均匀化理论和双线性内...
  • 使用边界值方法设计测试用例笔记

    千次阅读 2019-09-27 03:30:18
    案例1:两位整数加法计算器 1、边界值的应用场合 ...一般情况下,需要对边界值(-99和99)以及边界值两边的数(-100和-98以及100和98)分别进行测试。 2、如何使用 把边界值的点(3个点)单独写用例 案例...
  • 程序优化

    万次阅读 2016-01-17 18:29:48
    程序优化是指利用软件开发工具对程序进行调整和改进,让程序充分利用资源,提高运行效率,缩减代码尺寸的过程。按照优化的侧重点不同,程序优化可分为运行速度优化和代码尺寸优化。 运行速度优化是指在充分掌握软...
  • C#与.NET3.5高级程序设计第四版高清PDF中文完整版

    千次下载 热门讨论 2011-07-05 10:25:50
    3.8 条件结构和关系/相等运算符  3.9 小结  第4章 c#核心编程结构ⅱ  4.1 方法和参数修饰符  4.2 成员重载  4.3 c#中的数组操作  4.4 枚举类型  4.5 结构类型  4.6 值类型和引用类型  4.7 值类型...
  • 如果晶格兼容,请考虑周期性边界条件。 (可选)打印移动的距离和距离矢量。 vasputil_chgarith CHG / CHGCAR格式文件中电荷密度的简单算法。 支持+,-,*,/用于元素算术运算,支持'avg'用于计算平均值。 vasputil_...
  • 本书全面介绍了delphi程序开发所用到的技术和技巧,共分19章,内容包括窗体与界面设计、控件应用、数据处理技术、图形技术、多媒体技术、文件系统、操作系统与window相关程序、注册表、数据库技术、sql查询相关技术...
  • 从计算机程序出现的第一天起,对效率的追求就是程序天生的坚定信仰,这个过程犹如一场没有终点,永不停歇的F1方程式竞赛,程序员试车手,技术平台则是在赛道上飞驰的赛车。 文章目录早期(编译期)优化概述Javac...
  • 项目管理和软件开发的边界

    千次阅读 2017-10-21 16:21:21
    程序员的人生就是和一个个的软件项目打交道的人生。 不能管理好项目过程的程序员不是好的开发人员。...所以如何区分好项目管理和软件开发的边界,统筹好二者才是一个好的互联网项目管理和软件开发过程。
  • 【微信小程序】从入门到放弃

    万次阅读 2016-10-25 21:40:52
    关于微信小程序是什么,能做什么的问题,草民在此不在罗列了,随着小程序的天天刷屏,想必您也是来吃一些干货,本篇博文和大家走进微信小程序的从入门到放弃~
  • PostSharp的边界方法
  • 数组的FFTView使用周期性边界条件进行索引,并将数组的所有索引向下移动1。 用法 让我们创建一个随机信号: julia > using FFTViews julia > a = rand ( 8 ) 8 - element Array{Float64, 1 } : 0.720657 0.42337...
  • 测试用例--等价类划分、边界值法

    千次阅读 2020-09-15 15:27:30
     2)边界值法()  3)因果图法  4)判定表法  5)正交排列法  6)测试大纲法  7)场景法(*****)  至少要掌握每种方法的适用场合(用在哪)和使用步骤(怎么用)  编写测试用例可以参考
  • NaSch模型的matlab实现

    2015-10-27 09:57:05
    程序根据1992年Nagel和Schreckenberg发表的交通流论文,利用matlab编写。采用周期性边界条件,希望有所帮助
  • 程序将分割图像作为输入并生成有限元模型来求解通过材料的静态热流,可以指定周期性边界条件,输出是模型的热导率。
  • 是把程序的输入域划分成若干个子集合(等价类),然后从每个子集合(等价类)中选取少数具有代表的数据作为测试的输入数据。 在该子集合中,所有的输入数据对于揭露软件中的错误都是等效的。-----减少测试用例数量,...
  • 作为一名程序开发者,应该努力使渲染引擎的刷新率维持在60fps,也就是说在每帧之间大约有16ms,这段时间包括了基本图元在图形硬件上的描画。具体内容如下: >尽可能的使用异步事件驱动来编程。 >使用工作
  • 译《Office商业应用程序入门》

    千次阅读 2008-06-19 16:56:00
    第一章:Office商业应用程序入门-罗伯·巴克,微软公司概述在过去20年里,公司和组织已花费了数十亿美元购买,安装,部署和维护line-of-business(LOB)系统来管理客户资料,库存,帐单,产品的生命周期,和许多其他...
  • 程序根据1992年Nagel和Schreckenberg发表的交通流论文,利用matlab编写。采用周期性边界条件,答案在array矩阵中,希望有所帮助

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 34,197
精华内容 13,678
关键字:

周期性边界条件程序