Fipy 在求解时间步长之间更改 CellValue

如何解决Fipy 在求解时间步长之间更改 CellValue

我想将网格（偏微分方程）与复杂的混合反应器（常微分方程）连接起来。我想取一点值，求解 ODE，然后重新输入网格中的值（将在其中求解 PDE）。

Clang

但是，我看到 .setValue 之后的变量没有用于新的迭代。我知道由于扩散而存在一些差异，但我认为这是不正确的：

 # Geometry
Lx = 2   # meters
Ly = 2    # meters
nx = 101 #int((Lx / 0.05) + 1) # nodes
ny = 101 #int((Ly / 0.05) + 1)

# Build the mesh:
mesh = Grid2D(Lx = Lx,Ly = Ly,nx = nx,ny = ny)
X,Y = mesh.cellCenters
xx,yy = mesh.faceCenters

C0_DO = 10. #mg/L # dissolved Oxygen initial concentration
C0_OM = 5.
T0 = 28. #mg/L # # Temperature
print(f'DO saturation concentration: {calc_DOsat(T0)} mg/L')

# dissolved oxygen:
C_DO = CellVariable(name="concentration_DO",mesh=mesh,value=C0_DO,#hasOld=True
                 )

# Temperature
T = CellVariable(name="temperature",value=T0,#hasOld=True
                 )

C_OM = CellVariable(name="concentration_OM",value=C0_OM,#hasOld=True
                 )

# Transport parameters (diffusion constants)

# dissolved oxygen diffusion constant
D_DO = 2. #m2/s
D_DO = FaceVariable(name='DO_diff',value=D_DO) # This is necessary to set a fixed-flux boundary condition
D_DO.constrain(0.,mesh.exteriorFaces) # This is necessary to set a fixed-flux boundary condition


D_OM = 2.

# Reaction and source terms

O_R = 1. # Re-aeration temperature correction factor (-)
Y_R = 1. #reaeration_rate() # Re-aeration rate (1/day)
DOsat = (14.652 - 0.41022 * T + 0.007991 * T ** 2 - 0.000077774 * T ** 3) # Oxygen saturation concentration (mg/L)
xi_SOD = 0.1 # Sediment Oxygen Demand (BC) (1/day)

# Oxygen transport-reaction
eqDO = (TransientTerm(var = C_DO) == # Transient term
        DiffusionTerm(coeff = D_DO,var = C_DO) # Diffusion term
        + (mesh.facesTop * (Y_R * (O_R ** (T.faceValue - 20)) * ((14.652 - 0.41022 * T.faceValue + 0.007991 * T.faceValue ** 2 - 0.000077774 * T.faceValue ** 3) - C_DO.faceValue))).divergence #Top Boundary Condition (fixed-flux)
        + (mesh.facesBottom * (xi_SOD * C_DO.faceValue)).divergence # Bottom Boundary Condition (fixed-flux)
        )

eqOM = (TransientTerm(var = C_OM) == # Transient term
        DiffusionTerm(coeff = D_OM,var = C_OM) # Diffusion term
        )

eqTOT = eqDO & eqOM


    def reactor_model(C,t):
      k1 = 0.5
      k2 = 0.1
    
      C1 = C[0]
      C2 = C[1]
    
      dC1dt = -k1 * C1
      dC2dt = -k1 * C1
    
      dCdt = [dC1dt,dC2dt]
      
      return dCdt



  # PDESolver hyperparameters
steps = 10
sweeps = 2
dt = 1e-1
reactor_act = True

t = time.time()
for step in range(steps):
  C_DO.updateOld()
  C_OM.updateOld()
  
  for sweep in range(sweeps):
    eqTOT.sweep(dt=dt)

  if reactor_act:
    # Reactor
    r_in = [C_DO([(Lx/2),(0)]).mean(),# dissolved Oxygen reactor input
            C_OM([(Lx/2),(0)]).mean()  # Organic Matter reactor input
            ]

    print('Reactor input:',r_in)

    # solve ODE
    C = odeint(reactor_model,r_in,np.linspace(0,dt,10))
    r_out = C[-1]

    print('Reactor output:',r_out)

    # Reinput ODE
    C_DO.setValue(r_out[0],where=(Y >= Ly - mesh.dy / 2 ) & (X <= (Lx/2 + mesh.dx / 2)) & (X > (Lx/2 - mesh.dx / 2)))
    C_OM.setValue(r_out[1],where=(Y >= Ly - mesh.dy / 2 ) & (X <= (Lx/2 + mesh.dx / 2)) & (X > (Lx/2 - mesh.dx / 2)))
  
    C_DO.setValue(r_out[0],where=(Y <= mesh.dy / 2 ) & (X <= (Lx/2 + mesh.dx / 2)) & (X > (Lx/2 - mesh.dx / 2)))
    C_OM.setValue(r_out[1],where=(Y <= mesh.dy / 2 ) & (X <= (Lx/2 + mesh.dx / 2)) & (X > (Lx/2 - mesh.dx / 2)))

elapsed = time.time() - t
print('Calculation completed')
print(f'Runtime: {elapsed/60} minutes')

有什么方法可以在迭代之间更改 CellValues 吗？

谢谢。