问题描述
首先让我说我在 NARKIVE FiPy 邮件列表存档中发现了类似的问题，但由于方程式无法加载，因此它们不是很有用。例如 Convection-diffusion problem on a 1D cylindrical grid，或在另一个邮件列表存档 Re: FiPy Heat Transfer Solution。在第二个链接邮件中，丹尼尔说：
在 FiPy 中有两种方法可以在圆柱域上求解。你可以 使用笛卡尔坐标中的标准扩散方程（第二个方程 下面）和一个实际上是圆柱形的网格，或者你可以使用 在圆柱坐标系上制定的扩散方程（第一 下面的方程）并使用标准的 2D / 1D 网格。
方程也不存在。在这种情况下，它实际上很好，因为我理解第一个解决方案并且我想使用它。
我想在一维圆柱网格上求解以下方程（对不起，我还没有 10 声望，所以我不能发布漂亮的渲染方程）：
有边界条件：
其中 rho_core 是网格的左侧，rho_edge 是网格的右侧。 rho 是归一化半径，J 是雅可比：
R 是以米为单位的真实半径，因此雅可比矩阵的维数是距离。初始条件并不重要，但在我的代码示例中，我将在 R=0.8 处使用数字 Dirac-delta。
我有一个没有（！）雅可比行列式的工作示例，但它很长，而且它不使用 FiPy 的查看器，所以我将链接一个要点：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)
我的代码可以使用与上述相同的边界条件求解以下方程：
为了简单起见，我使用空间独立系数，唯一的例外是在内边缘，其中对流系数为 0，扩散系数几乎为 0。在链接的代码中，我使用了均匀分布的初始条件。
我的第一个问题是，为什么在使用 fipy.Grid1D 和 fipy.CylindricalGrid1D 时会得到完全相同的结果？我应该得到不同的结果，对吧？我应该如何重写我的代码才能区分简单的一维网格和一维圆柱网格？
我的实际问题不在于这个确切的代码，我只是想简化我的问题，但正如评论中所指出的，这段代码在不同的网格中不会产生相同的结果。所以我只会发布一个指向 Jupyter Notebook 的 GitHub 链接，它在未来可能会停止工作。
The Jupyter Notebook 如果你想运行它，第一个代码单元应该首先运行，然后只有最后一个单元是相关的。忽略参考图像。线图显示了扩散和对流系数。当我使用 Grid1D 或 CylindricalGrid1D 运行最后一个单元格时，我得到了相同的结果（我非常精确地比较了这些图）
抱歉，我无法重命名所有变量，所以我希望根据我的评论和上面更改的代码（我也更改了代码中的注释），您可以理解我要做什么。
>我的另一个问题是关于雅可比矩阵。我该如何实施？我查看了文档中使用雅可比行列式的唯一示例，但雅可比行列式是一个矩阵，并且还使用了 scipy.optimize.fsolve() 函数。
解决方法
[从评论中的讨论中拼凑出一个答案]
Grid1D 和 CylindricalGrid1D 之间的结果相似，尤其是在早期步骤中，但它们并不相同。随着问题的发展，它们是完全不同的。
FiPy 不喜欢发散之外的东西，但您应该能够将方程乘以 J 并将其放入 TransientTerm 的系数中，例如，
或
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 返回一个向量，即使是一维的，但系数必须是标量，所以取点积得到 x 分量。