运行时错误:Factor完全是单数

Question

我正在尝试实现2温度模型

1. C_e(∂T_e)/∂t=∇[k_e∇T_e ]-G(T_e-T_ph )+ A(r,t)

2. C_ph(∂T_ph)/∂t=∇[k_ph∇T_ph] + G(T_e-T_ph)

代码

from fipy.tools import numerix
import  scipy
import fipy
import numpy as np
from fipy import CylindricalGrid1D
from fipy import Variable, CellVariable, TransientTerm, DiffusionTerm, Viewer, LinearLUSolver, LinearPCGSolver, 
LinearGMRESSolver, ImplicitDiffusionTerm, Grid1D
FIPY_SOLVERS = scipy
## Mesh

nr = 50
dr = 1e-7
# r = nr * dr
mesh = CylindricalGrid1D(nr=nr, dr=dr, origin=0)
x = mesh.cellCenters[0]
# Variables
T_e = CellVariable(name="electronTemp", mesh=mesh,hasOld=True)
T_e.setValue(300)
T_ph = CellVariable(name="phononTemp", mesh=mesh, hasOld=True)
T_ph.setValue(300)
G = CellVariable(name="EPC", mesh=mesh)
t = Variable()
# Material parameters
C_e = CellVariable(name="C_e", mesh=mesh)
k_e = CellVariable(name="k_e", mesh=mesh)
C_ph = CellVariable(name="C_ph", mesh=mesh)
k_ph = CellVariable(name="k_ph", mesh=mesh)
C_e = 4.15303 - (4.06897 * numerix.exp(T_e / -85120.8644))
C_ph = 4.10446 - 3.886 * numerix.exp(-T_ph / 373.8)
k_e = 0.1549 * T_e**-0.052
k_ph =1.24 + 16.29 * numerix.exp(-T_ph / 151.57)

G = numerix.exp(21.87 + 10.062 * numerix.log(numerix.log(T_e )- 5.4))
# Boundary conditions
T_e.constrain(300, where=x > 4.5e-6)
T_ph.constrain(300, where=x > 4.5e-6)
# Source  𝐴(𝑟,𝑡) = 𝑎𝐷(𝑟)𝜏−1 𝑒−𝑡/𝜏   , 𝐷(𝑟) = 𝑆𝑒 exp (−𝑟2/𝜎2)/√2𝜋𝜎2
sig = 1.0e-6
tau = 1e-15
S_e = 35
d_r = (S_e * 1.6e-9 * numerix.exp(-x**2 /sig**2)) / (numerix.sqrt(2. * 3.14 * sig**2))
A_t = numerix.exp(-t/tau)
a = (numerix.sqrt(2. * 3.14)) / (3.14 * sig)
A_r = a * d_r * tau**-1 * A_t 

eq0 = (TransientTerm(var=T_e, coeff=C_e) == DiffusionTerm(var=T_e, coeff=k_e) - G*(T_e - T_ph) + A_r
eq1 =(TransientTerm(var=T_ph, coeff=C_ph) == DiffusionTerm(var=T_ph, coeff=k_ph) + G*(T_e - T_ph)
eq = eq0 & eq1
dt = 1e-18
steps = 7000
elapsed = 0.
vi = Viewer((T_e, T_ph), datamin=0., datamax=2e4)

for step in range(steps):
T_e.updateOld()
T_ph.updateOld()
vi.plot()
res = 1e100
dt *= 1.1
while res > 1:
res = eq.sweep(dt=dt)
print(t, res)
t.setValue(t + dt)

问题代码工作很好，非常小的dt = 1e-18，但我需要运行它，直到e 1e-10。

使用这个时间步骤将花费很长时间，并且设置dt *= 1.1，结果在某个点开始增加，然后

给出以下运行时错误:

因子是完全奇异的

即使是非常小的增量dt*= 1.005也会出现同样的问题。使用dt= 1.001运行代码很长一段时间，然后剩余会卡在某个值上。

问题

1. 这些方程的形式是否有错误?
2. 错误的原因是什么?
3. 是由于时间步长增加导致的错误吗?如果是，如何增加我的时间步长?

Answer 1

我对代码做了一些更改，可以使您的运行时间为1e-10。主要变化如下

• 对于G的条款使用ImplicitSourceTerm。这使溶液稳定。
• 在扫描步骤中应用underRelaxation=0.5。这减慢了扫描循环中的更新，因此反馈循环被抑制了。
删除FIPY_SOLVERS=scipy。这没有做任何事情。FIPY_SOLVERS是一个环境变量，可以在Python环境之外设置。
• 应用边界条件的方式似乎很奇怪，所以我以更规范的方式应用它们。
• 扫描循环固定为10次扫描，以快速达到稳定状态。注意，当解接近稳定状态时，残差不一定会变好。如果你需要一个准确的瞬态，可能需要回到剩余检查。

from fipy.tools import numerix
import  scipy
import fipy
import numpy as np
from fipy import CylindricalGrid1D
from fipy import Variable, CellVariable, TransientTerm, DiffusionTerm, Viewer, LinearLUSolver, LinearPCGSolver, 
LinearGMRESSolver, ImplicitDiffusionTerm, Grid1D, ImplicitSourceTerm
## Mesh

nr = 50
dr = 1e-7
# r = nr * dr
mesh = CylindricalGrid1D(nr=nr, dr=dr, origin=0)
x = mesh.cellCenters[0]
# Variables
T_e = CellVariable(name="electronTemp", mesh=mesh,hasOld=True)
T_e.setValue(300)
T_ph = CellVariable(name="phononTemp", mesh=mesh, hasOld=True)
T_ph.setValue(300)
G = CellVariable(name="EPC", mesh=mesh)
t = Variable()
# Material parameters
C_e = CellVariable(name="C_e", mesh=mesh)
k_e = CellVariable(name="k_e", mesh=mesh)
C_ph = CellVariable(name="C_ph", mesh=mesh)
k_ph = CellVariable(name="k_ph", mesh=mesh)
C_e = 4.15303 - (4.06897 * numerix.exp(T_e / -85120.8644))
C_ph = 4.10446 - 3.886 * numerix.exp(-T_ph / 373.8)
k_e = 0.1549 * T_e**-0.052
k_ph =1.24 + 16.29 * numerix.exp(-T_ph / 151.57)

G = numerix.exp(21.87 + 10.062 * numerix.log(numerix.log(T_e )- 5.4))
# Boundary conditions
T_e.constrain(300, where=mesh.facesRight)
T_ph.constrain(300, where=mesh.facesRight)
# Source  𝐴(𝑟,𝑡) = 𝑎𝐷(𝑟)𝜏−1 𝑒−𝑡/𝜏   , 𝐷(𝑟) = 𝑆𝑒 exp (−𝑟2/𝜎2)/√2𝜋𝜎2
sig = 1.0e-6
tau = 1e-15
S_e = 35
d_r = (S_e * 1.6e-9 * numerix.exp(-x**2 /sig**2)) / (numerix.sqrt(2. * 3.14 * sig**2))
A_t = numerix.exp(-t/tau)
a = (numerix.sqrt(2. * 3.14)) / (3.14 * sig)
A_r = a * d_r * tau**-1 * A_t

eq0 = (
TransientTerm(var=T_e, coeff=C_e) == 
DiffusionTerm(var=T_e, coeff=k_e) - 
ImplicitSourceTerm(coeff=G, var=T_e) + 
ImplicitSourceTerm(var=T_ph, coeff=G) + 
A_r)
eq1 = (TransientTerm(var=T_ph, coeff=C_ph) == DiffusionTerm(var=T_ph, coeff=k_ph) + ImplicitSourceTerm(var=T_e, coeff=G) - ImplicitSourceTerm(coeff=G, var=T_ph))
eq = eq0 & eq1
dt = 1e-18
steps = 7000
elapsed = 0.
vi = Viewer((T_e, T_ph), datamin=0., datamax=2e4)

for step in range(steps):
T_e.updateOld()
T_ph.updateOld()
vi.plot()
res = 1e100
dt *= 1.1
count = 0
while count < 10:
res = eq.sweep(dt=dt, underRelaxation=0.5)
print(t, res)
count += 1
print('elapsed:', t.value)
t.setValue(t + dt)

关于你的问题。

这些方程的形式是否有错误?

事实上,没有。形式主义没有错，但最好使用ImplicitSourceTerm。

是什么导致错误?

在这个系统中有两个不稳定的来源。当明确地写出来时，方程中的源项在某个时间步长以上是不稳定的。使用ImplcitSourceTerm可以消除这种不稳定性。在方程的耦合中也存在某种不稳定性。

错误是因为时间步长增加吗?如果是，如何增加我的时间步长?

上面解释道。

Answer 2

除了@wd15的回答:

你的方程是非常非线性的。您可能会受益于牛顿迭代，以获得良好的收敛性。

正如@TimRoberts所说，不加限制地几何增加时间步长可能不是一个好主意。

我最近发布了一个名为steppyngstones的包，它负责调整时间步长。尽管它是一个独立的包，但它打算与FiPy一起使用。例如，您可以将solve循环更改为:

from steppyngstounes import FixedStepper, PIDStepper
T_e.updateOld()
T_ph.updateOld()
for checkpoint in FixedStepper(start=0, stop=1e-10, size=1e-12):
for step in PIDStepper(start=checkpoint.begin,
stop=checkpoint.end,
size=dt):
res = 1e100
for sweep in range(10):
res = eq.sweep(dt=dt, underRelaxation=0.5)
print(t, sweep, res)
if step.succeeded(error=res / 1000):
T_e.updateOld()
T_ph.updateOld()
t.value = step.end
else:
T_e.value = T_e.old
T_ph.value = T_ph.old
print('elapsed:', t.value)
# the last step might have been smaller than possible,
# if it was near the end of the checkpoint range
dt = step.want
_ = checkpoint.succeeded()
vi.plot()

此代码将每隔1e-12个时间单位更新查看器，并自适应地在这些检查点之间进行切换。如果让事情变得更有趣的话，包中还有其他的步进器可以方便地以几何或指数方式增加检查点。

通过减少扫描次数并让适配器在开始时采取更小的时间步，您可能会获得更好的整体性能。我发现没有一个时间步长小到足以使初始残差低于777.9。在最初的几个步骤之后，误差度量可能会更加激进，给出更准确的结果。