pySPH学习

NightFuryFFF

已于 2024-08-01 20:40:39 修改

阅读量701

点赞数 7

文章标签：学习 python

于 2024-07-24 10:00:04 首次发布

本文链接：https://blog.csdn.net/NightFuryFFF/article/details/140549435

版权

example taylor_bar 2D

主程序代码：mian.py

"""Taylor bar example with SPH. (5 minutes)
"""

# 导入函数
import numpy

from pysph.sph.equation import Group

from pysph.base.utils import get_particle_array
from pysph.base.kernels import WendlandQuintic
from pysph.solver.application import Application
from pysph.solver.solver import Solver
from pysph.sph.integrator import PECIntegrator
from pysph.sph.integrator_step import SolidMechStep

# basic sph equations
from pysph.sph.basic_equations import ContinuityEquation, \
    MonaghanArtificialViscosity, XSPHCorrection, VelocityGradient2D

# baic stress equations
from pysph.sph.solid_mech.basic import HookesDeviatoricStressRate,\
    MomentumEquationWithStress, EnergyEquationWithStress

# plasticity model and eos
from pysph.sph.solid_mech.hvi import VonMisesPlasticity2D, MieGruneisenEOS

# boundary force
from pysph.sph.boundary_equations import MonaghanBoundaryForce

# Numerical Parameters and constants
# 数值参数和常数
dx = dy = 0.000384848
hdx = 2.0
h = hdx * dx
r0 = 7850.0
m0 = dx * dy * r0
v_s = 200.0
ss = 4699.0
C = 3630.0
S = 1800.0
gamma = 1.81
alpha = 0.5
beta = 0.5
eta = 0.01
eps = 0.5
bar_width = 0.0076
G = 8*1e10
Yo = 6*1e8
ro2 = 2750.0
plate_start = -2.0*bar_width
plate_end = 2.0*bar_width

# 粒子初始化函数

# 创建平板粒子
# 初始化粒子属性，如法线、切线、声速和质量
def get_plate_particles():
    x = numpy.arange(plate_start, plate_end+dx, dx)
    y = numpy.zeros_like(x)

    # normals and tangents
    tx = numpy.ones_like(x)
    ty = numpy.zeros_like(x)
    tz = numpy.zeros_like(x)

    ny = numpy.ones_like(x)
    nx = numpy.zeros_like(x)
    nz = numpy.zeros_like(x)
    cs = numpy.ones_like(x)*ss
    pa = get_particle_array(name='plate', x=x, y=y, tx=tx, ty=ty, tz=tz,
                            nx=nx, ny=ny, nz=nz, cs=cs)
    pa.m[:] = m0

    return pa

# 创建矩形杆粒子
# 初始化粒子属性，如平滑长度、质量、密度、声速和初始速度
def get_bar_particles():
    xarr = numpy.arange(-bar_width/2.0, bar_width/2.0 + dx, dx)
    yarr = numpy.arange(4*dx, 0.0254 + 4*dx, dx)

    x,y = numpy.meshgrid( xarr, yarr )
    x, y = x.ravel(), y.ravel()

    print('Number of bar particles: ', len(x))

    hf = numpy.ones_like(x) * h
    mf = numpy.ones_like(x) * dx * dy * r0
    rhof = numpy.ones_like(x) * r0
    csf = numpy.ones_like(x) * ss
    z = numpy.zeros_like(x)
    pa = get_particle_array(name="bar",
                            x=x, y=y, h=hf, m=mf, rho=rhof, cs=csf,
                            e=z)
    # negative fluid particles
    pa.v[:]=-200
    return pa

class TaylorBar(Application):
    def create_particles(self):
        bar = get_bar_particles()
        plate = get_plate_particles()

        # add requisite properties

        # velocity gradient for the bar
        for name in ('v00', 'v01', 'v02', 'v10', 'v11', 'v12', 'v20', 'v21', 'v22'):
            bar.add_property(name)

        # deviatoric stress components
        for name in ('s00', 's01', 's02', 's11', 's12', 's22'):
            bar.add_property(name)

        # deviatoric stress accelerations
        for name in ('as00', 'as01', 'as02', 'as11', 'as12', 'as22'):
            bar.add_property(name)

        # deviatoric stress initial values
        for name in ('s000', 's010', 's020', 's110', 's120', 's220'):
            bar.add_property(name)

        bar.add_property('e0')

        # artificial stress properties
        for name in ('r00', 'r01', 'r02', 'r11', 'r12', 'r22'):
            bar.add_property(name)

        # standard acceleration variables and initial values.
        for name in ('arho', 'au', 'av', 'aw', 'ax', 'ay', 'az', 'ae',
                     'rho0', 'u0', 'v0', 'w0', 'x0', 'y0', 'z0', 'e0'):
            bar.add_property(name)

        bar.add_constant('G', G)
        bar.add_constant('n', 4)

        kernel = WendlandQuintic(dim=2)
        wdeltap = kernel.kernel(rij=dx, h=hdx*dx) #核函数的值，为常量
        bar.add_constant('wdeltap', wdeltap)

        return [bar, plate]

    def create_solver(self):
        kernel = WendlandQuintic(dim=2)
        self.wdeltap = kernel.kernel(rij=dx, h=hdx*dx) # 核函数的值，为常数

        integrator = PECIntegrator(bar=SolidMechStep())

        solver = Solver(kernel=kernel, dim=2, integrator=integrator)
        dt = 1e-9
        tf = 2.5e-5
        solver.set_time_step(dt)
        solver.set_final_time(tf)
        return solver

    def create_equations(self):
        equations = [

            # Properties computed set from the current state
            Group(
                equations=[
                    # p
                    MieGruneisenEOS(dest='bar', sources=None, gamma=gamma,
                                    r0=r0, c0=C, S=S),

                    # vi,j : requires properties v00, v01, v10, v11
                    VelocityGradient2D(dest='bar', sources=['bar',]),

                    # rij : requires properties s00, s01, s11
                    VonMisesPlasticity2D(flow_stress=Yo, dest='bar',
                                         sources=None),
                    ],
                ),

            # Acceleration variables are now computed
            Group(
                equations=[

                    # arho
                    ContinuityEquation(dest='bar', sources=['bar']),

                    # au, av
                    MomentumEquationWithStress(
                        dest='bar', sources=['bar']),

                    # au, av
                    MonaghanArtificialViscosity(
                        dest='bar', sources=['bar'], alpha=0.5, beta=0.5),

                    # au av
                    MonaghanBoundaryForce(
                        dest='bar', sources=['plate'], deltap=dx),

                    # ae
                    EnergyEquationWithStress(dest='bar', sources=['bar'],
                                             alpha=0.5, beta=0.5, eta=0.01),

                    # a_s00, a_s01, a_s11
                    HookesDeviatoricStressRate(
                        dest='bar', sources=None),

                    # ax, ay, az
                    XSPHCorrection(
                        dest='bar', sources=['bar',], eps=0.5),

                    ]

                ) # End Acceleration Group

        ] # End Group list
        return equations


if __name__ == '__main__':
    app = TaylorBar()
    app.run()

总体分析

时间积分框架示意

PySPH参考文档

此处翻译并记录

PySPH implementation

example dam_break.2d

two ParticleArrays 两组粒子数据

pysph.base.nnps.NNPS 邻域搜索

Specifying the equations 指定方程

给定粒子数组，我们通过向求解器(见pysph.solver.solver.Solver)传递一个Equation对象列表(见SPH方程)来请求对粒子执行一组给定的操作。

equations = [
    # Equation of state
    Group(equations=[
        ...
    ], ),

    # Acceleration
    Group(equations=[
        ...
    ] ),
]

python笔记：Class及其实例化

python类的定义与实例化-CSDN博客

我们看到我们使用了两个Group对象(参见pysph.sph.equation.Group)，分离了逻辑依赖的求值的两个部分。计算加速度的第二组必须在更新压力的第一组之后进行计算。

class pysph.sph.equation.Group

classpysph.sph.equation.Group(equations, real=True, update_nnps=False, iterate=False, max_iterations=1, min_iterations=0, pre=None, post=None, condition=None, start_idx=0, stop_idx=None, name=None)

source:

real (bool) - 指定是否只对非远程/非鬼影粒子进行操作。

Notes

当并行运行模拟时，通常应该在每个处理器中运行所有粒子(本地和远程)的总和密度。这是因为我们必须在当前处理器中更新远程邻居的压力/密度。否则，由于远端粒子具有较老的密度，结果可能不正确。TaitEOS也是如此。在这些情况下，计算方程的组应该将real设为False。

回想一下，在我们假设的实现中，我们必须对SPHCalc.compute方法进行类似的隔离:

class SPHCalc:
    def __init__(nnps, particles):
        ...

    def compute(self):
        self.eos()
        self.accelerations()

PySPH将尊重用户提供的方程和方程组的顺序。这种灵活性也意味着很容易犯细微的错误。

一般来说，SPH算法的过程如下伪代码所示:

for destination in particles:
    for equation in equations:
        equation.initialize(destination)

# This is where bulk of the computation happens.
for destination in particles:
    for source in destination.neighbors:
        for equation in equations:
            equation.loop(source, destination)

for destination in particles:
    for equation in equations:
        equation.post_loop(destination)

# Reduce any properties if needed.
total_mass = reduce_array(particles.m, 'sum')
max_u = reduce_array(particles.u, 'max')

对每一个目标粒子，循环每一个方程，初始化方程

对每一个目标粒子，循环每一个源粒子，再循环计算每一个方程的loop

使用最近邻算法识别给定粒子的邻居。PySPH会自动为用户完成这项工作，并在内部使用基于链表的算法来识别邻居。

在PySPH中，我们在写方程时遵循一些简单的约定。让我们先看几个方程。为了与我们假设的实现和SPHCalc保持类比，我们考虑了 pysph.sph.basic.TaitEOS 和 pysph.sph.basic_ContinuityEquation 对象的实现。

前者看起来像:

MomentumEquation

细看代码结构

if __name__ == '__main__':
    app = TaylorBar()
    app.run()

run()

run() 方法是应用程序的主要入口点，负责初始化和运行主要功能，调用一次
包含 setup() 和 solve() ，分别进行设置和求解。三者均在 application.py 中
argv 参数使得应用程序可以在命令行环境和交互环境下灵活运行。

solve()

求解 in application.py，class Application(object):

def solve(self):
        """This runs the solver.

        Note that this method expects that `setup has already been called.

        Don't use this method unless you really know what you are doing.
        """
        start_time = time.time()

        # 运行求解器
        self.solver.solve(not self.options.quiet)

        end_time = time.time()
        run_duration = end_time - start_time
        self._message("Run took: %.5f secs" % (run_duration))
        self._write_info(
            self.info_filename, completed=True, cpu_time=run_duration
        )

        self._stop_interfaces()
        self._write_profile_info()

运行求解器就一行，循环在这其中进行：

self.solver.solve(not self.options.quiet)

1. self.solver 对象
定义：self.solver 是 Application 类的一个实例变量，代表一个求解器对象。
职责：求解器对象负责执行模拟计算的实际过程。它可能是 pysph.solver.solver.Solver 类的一个实例，或者是某个自定义求解器类的实例。
作用：该求解器对象封装了模拟计算中涉及的所有数值方法、迭代过程和求解逻辑。
2. solve 方法
定义：solve 是 self.solver 对象中的一个方法。
职责：solve 方法启动求解器的计算过程，进行流体动力学模拟。
作用：该方法负责按设定的方案和初始条件推进模拟至完成。
3. 参数 not self.options.quiet
定义：self.options.quiet 是 Application 类中的一个布尔属性，通常从命令行参数或默认设置中获取。

职责：该属性控制输出日志的详细程度。

作用：

self.options.quiet 为 True：求解器在运行时不输出详细信息。
self.options.quiet 为 False：求解器在运行时输出详细信息，便于调试和监控进程。
参数传递：通过传递 not self.options.quiet，代码控制日志输出的开关。

self.solver

in setup(), self.solver = self.create_solver()

在主程序中：

    def create_solver(self):
        # 设置核函数
        kernel = WendlandQuintic(dim=2)
        self.wdeltap = kernel.kernel(rij=dx, h=hdx*dx)

        # 积分方法
        integrator = PECIntegrator(bar=SolidMechStep())

        # 导入求解器
        solver = Solver(kernel=kernel, dim=2, integrator=integrator)
        # 设置时间步长、总时间
        dt = 1e-9
        tf = 2.5e-5
        solver.set_time_step(dt)
        solver.set_final_time(tf)
        return solver

这个求解器是怎么样的？

Solver 的输入：核函数自由度时间积分方法：PECIntegrator(bar=SolidMechStep)

pysph/solver/solver.py class Solver(object)

python笔记：@函数装饰器 Python @函数装饰器及用法（超级详细） (biancheng.net)

使用函数装饰器 A() 去装饰另一个函数 B()，其底层执行了如下 2 步操作：
1 将 B 作为参数传给 A() 函数；
2 将 A() 函数执行完成的返回值反馈回 B()。

#funA 作为装饰器函数
def funA(fn):
    print("AAA")
    fn() # 执行传入的fn参数
    print("BBB")
    return "CCC"

@funA
def funB():
    print("DDD")


print(funB)

程序执行流程：

AAA
DDD
BBB

CCC

class SolidMechStep(IntegratorStep):

initialize

stage1

stage2

粒子属性存在名为 bar 的变量中，由 get_particle_array(additional_props=None, constants=None, backend=None, props): 设置

默认属性 + 额外添加属性

create_equations

在当前步内计算的参数

 # Properties computed set from the current state
            Group(
                equations=[
                    # p
                    MieGruneisenEOS(dest='bar', sources=None, gamma=gamma,
                                    r0=r0, c0=C, S=S),

                    # vi,j : requires properties v00, v01, v10, v11
                    VelocityGradient2D(dest='bar', sources=['bar',]),

                    # rij : requires properties s00, s01, s11
                    VonMisesPlasticity2D(flow_stress=Yo, dest='bar',
                                         sources=None),
                    ],
                ),

1 MieGruneisenEOS

状态方程，根据密度、内能计算压力

2 VelocityGradient2D

计算速度梯度

输入：目标粒子及源粒子的位置、速度分量

输出：目标粒子的速度梯度

3 VonMisesPlasticity2D

屈服准则，屈服强度Y是直接输入的 Yo

检查材料的应力状态是否超过了指定的屈服应力，如果超过了，则进行塑性修正以确保应力不会超过屈服极限。

输入：直接使用dest的数据并更新，一个粒子索引及其应力张量的不同分量

输出：更新后的应力张量

计算加速度的公式组

计算加速度的第二组必须在更新压力的第一组之后进行计算。

# Acceleration variables are now computed
            Group(
                equations=[

                    # arho
                    ContinuityEquation(dest='bar', sources=['bar']),

                    # au, av
                    MomentumEquationWithStress(
                        dest='bar', sources=['bar']),

                    # au, av
                    MonaghanArtificialViscosity(
                        dest='bar', sources=['bar'], alpha=0.5, beta=0.5),

                    # au av
                    MonaghanBoundaryForce(
                        dest='bar', sources=['plate'], deltap=dx),

                    # ae
                    EnergyEquationWithStress(dest='bar', sources=['bar'],
                                             alpha=0.5, beta=0.5, eta=0.01),

                    # a_s00, a_s01, a_s11
                    HookesDeviatoricStressRate(
                        dest='bar', sources=None),

                    # ax, ay, az
                    XSPHCorrection(
                        dest='bar', sources=['bar',], eps=0.5),

                    ]

                ) # End Acceleration Group

1 ContinuityEquation

class ContinuityEquation(Equation):
    r"""Density rate:

    :math:`\frac{d\rho_a}{dt} = \sum_b m_b \boldsymbol{v}_{ab}\cdot
    \nabla_a W_{ab}`

    """
    def initialize(self, d_idx, d_arho):
        d_arho[d_idx] = 0.0


    def loop(self, d_idx, d_arho, s_idx, s_m, DWIJ, VIJ):
        vijdotdwij = DWIJ[0]*VIJ[0] + DWIJ[1]*VIJ[1] + DWIJ[2]*VIJ[2]
        d_arho[d_idx] += s_m[s_idx]*vijdotdwij

初始化 + 循环

输入：目标粒子位置、源粒子位置及质量、粒子对的核函数梯度及速度差

输出：该粒子的速度对时间的导数