使用 FiPy 在一维圆柱坐标系上求解 PDE

Question

首先让我说一下，我在 NARKIVE FiPy 邮件列表存档中发现了与我类似的问题，但由于方程式无法加载，因此它们不是很有用。例如 Convection-diffusion problem on a 1D cylindrical grid, or on another mailing list archive Re: FiPy Heat Transfer Solution。在第二封 linked 邮件中丹尼尔说：

There are two ways to solve on a cylindrical domain in FiPy. You can either use the standard diffusion equation in Cartesian coordinates (2nd equation below) and with a mesh that is actually cylindrical in shape or you can use the diffusion equation formulated on a cylindrical coordinate system (1st equation below) and use a standard 2D / 1D grid mesh.

方程式不存在。在这种情况下它实际上很好，因为我理解第一个解决方案并且我想使用它。

我想在一维圆柱网格上求解以下方程（抱歉，我还没有 10 个声望，所以我不能 post 漂亮的渲染方程）：

$\frac{\partial n(\rho,t)}{\partial t}=- \frac{1}{J}\frac{\partial}{\partial \rho} \Big[JK(\rho,t)n(\rho,t)\Big]+\frac{1}{J}\frac{\partial^2}{\partial \rho^2} \Big[JD(\rho,t)n(\rho,t)\Big]$

有边界条件：

$\frac{\partial n(t)}{\partial \rho}\Bigg|_{\rho=\rho{core}}=0 \quad\text{and}\quad n\bigg|_{\rho=\rho{edge}}=0$

其中 rho_core 是网格的左侧，rho_edge 是网格的右侧。 rho 是归一化半径，J 是雅可比行列式：

$J=\frac{\partial R}{\partial \rho}$

R是以米为单位的真实半径，所以雅可比矩阵的量纲是距离。初始条件并不重要，但在我的代码示例中，我将在 R=0.8.

处使用数字 Dirac-delta

我有一个没有（！）雅可比行列式的工作示例，但它很长，而且它没有使用 FiPy 的查看器所以我将 link 一个要点：https://gist.github.com/leferi99/142b90bb686cdf5116ef5aee425a4736

问题主要部分如下：

    import fipy as fp  ## finite volume PDE solver
    from fipy.tools import numerix  ## requirement for FiPy, in practice same as numpy
    import copy  ## we need the deepcopy() function because some FiPy objects are mutable
    import numpy as np
    import math

    ## numeric implementation of Dirac delta function
    def delta_func(x, epsilon, coeff):
        return ((x < epsilon) & (x > -epsilon)) * \
            (coeff * (1 + numerix.cos(numerix.pi * x / epsilon)) / (2 * epsilon))

    rho_from = 0.7  ## normalized inner radius
    rho_to = 1.  ## normalized outer radius
    nr = 1000  ## number of mesh cells
    dr = (rho_to - rho_from) / nr  ## normalized distance between the centers of the mesh cells
    duration = 0.001  ## length of examined time evolution in seconds
    nt = 1000  ## number of timesteps
    dt = duration / nt  ## length of one timestep
    
    ## 3D array for storing the density with the correspondant normalized radius values
    ## the density values corresponding to the n-th timestep will be in the n-th line 
    solution = np.zeros((nt,nr,2))
    ## loading the normalized radial coordinates into the array
    for j in range(nr):
        solution[:,j,0] = (j * dr) + (dr / 2) + rho_from
    
    mesh = fp.CylindricalGrid1D(dx=dr, nx=nr)  ## 1D mesh based on the normalized radial coordinates 
    mesh = mesh + (0.7,)  ## translation of the mesh to rho=0.7
    n = fp.CellVariable(mesh=mesh)  ## fipy.CellVariable for the density solution in each timestep
    
    diracLoc = 0.8  ## location of the middle of the Dirac delta
    diracCoeff = 1.  ## Dirac delta coefficient ("height")
    diracPercentage = 2  ## width of Dirac delta (full width from 0 to 0) in percentage of full examined radius
    diracWidth = int((nr / 100) * diracPercentage)
    
    ## diffusion coefficient
    diffCoeff = fp.CellVariable(mesh=mesh, value=100.)
    ## convection coefficient - must be a vector
    convCoeff = fp.CellVariable(mesh=mesh, value=(1000.,))
    
    ## applying initial condition - uniform density distribution
    n.setValue(1)

    ## boundary conditions
    gradLeft = (0.,)  ## density gradient (at the "left side of the radius") - must be a vector
    valueRight = 0.  ## density value (at the "right end of the radius")
    n.faceGrad.constrain(gradLeft, where=mesh.facesLeft)  ## applying Neumann boundary condition
    n.constrain(valueRight, mesh.facesRight)  ## applying Dirichlet boundary condition
    convCoeff.setValue(0, where=mesh.x<(R_from + dr))  ## convection coefficient 0 at the inner edge
    
    ## the PDE
    eq = (fp.TransientTerm() == fp.DiffusionTerm(coeff=diffCoeff)
          - fp.ConvectionTerm(coeff=convCoeff))
    
    ## Solving the PDE and storing the data
    for i in range(nt):
        eq.solve(var=n, dt=dt)
        solution[i,0:nr,1]=copy.deepcopy(n.value)

我的代码可以使用与上述相同的边界条件求解以下方程：

$\frac{\partial n(\rho,t)}{\partial t}=- \frac{\partial}{\partial \rho} \Big[K(\rho)n(\rho,t)\Big]+\frac{\partial^2}{\partial \rho^2} \Big[D(\rho)n(\rho,t)\Big]$

为了简单起见，我使用了空间无关的系数，唯一的例外是在内边缘，对流系数为 0，扩散系数几乎为 0。在 linked 代码中，我使用均匀分布的初始条件。

我的第一个问题是为什么在使用 fipy.Grid1D 和 fipy.CylindricalGrid1D 时得到完全相同的结果？我应该得到不同的结果，对吧？我应该如何重写我的代码才能区分简单的一维网格和一维圆柱网格？

我的实际问题不在于这个确切的代码，我只是想简化我的问题，但正如评论中所指出的那样，这段代码不会对不同的网格产生相同的结果。所以我只会 post 一个 GitHub link 到一个 Jupyter Notebook，它可能会在未来停止工作。 The Jupyter Notebook 如果你想运行它，第一个代码单元格应该是运行，然后只有最后一个单元格是相关的。忽略参考图像。线图显示扩散和对流系数。当我运行最后一个带有 Grid1D 或 CylindricalGrid1D 的单元格时，我得到了相同的结果（我非常精确地比较了这些图）

抱歉，我无法重命名我的所有变量，所以我希望根据我的评论和上面更改的代码（我也更改了代码中的注释）您可以理解我正在尝试做什么。

我的另一个问题是关于雅可比矩阵的。我该如何实施？我查看了文档中唯一使用雅可比行列式的示例，但该雅可比行列式是一个矩阵，并且它还使用了 scipy.optimize.fsolve() 函数。

Answer 1

[根据评论中的讨论拼凑答案]

a Grid1D 和 a CylindricalGrid1D 的结果相似，尤其是在早期步骤中，但它们并不相同。随着问题的发展，它们是完全不同的。

FiPy 不喜欢散度以外的东西，但你应该可以将方程乘以 J 并将其放入 TransientTerm 的系数中，例如， $\frac{\partial J n(\rho,t)}{\partial t}=\nabla\cdot\left[J \vec{K}(\rho, t) n(\rho, t) \right ] + \nabla\cdot\left[J D(\rho, t) \nabla n(\rho, t)\right ]$

或

eq = fp.TransientTerm(J) == fp.DiffusionTerm(coeff=J * diffCoeff) - fp.ConvectionTerm(coef=J * convCoeff)

对于雅可比矩阵，您可以根据归一化半径为实际半径创建 CellVariable，然后取其梯度：

real_radius = fp.CellVariable(mesh=mesh, value=...)
J = real_radius.grad.dot([[1]])

.grad returns一个向量，即使是一维的，但系数必须是标量，所以取点积得到x分量。

使用 FiPy 在一维圆柱坐标系上求解 PDE

Solving PDE on 1D cylindrical coordinates with FiPy

python

fipy