I am solving the same differential equation using two different methods. One uses finite differences (code A) and the other uses a neural network in Nvidia Modulus, a deep learning framework built on PyTorch (code B). The problem is, that the differential equation in code A has a k term in def rhs , in the differential equation, and I do not know how to code for it in code B. Please help me in this area.
Below is my executable code A :
import numpy as np
nz = 50
dz = 1 / nz
nt = 50
t = np.linspace(0, 4, nt)
p = np.array([0.01, 0.5, 1.0])
#Initialize concentration and temperature arrays
C = np.zeros(nz)
T = np.zeros(nz)
#Set initial condition for concentration
C[0] = 1.0
#Initialize arrays to store results
c_results = np.zeros((nt, nz))
T_results = np.zeros((nt, nz))
#Define the Arrhenius rate constant function
def k(temperature):
arrh_expr = 4.68 * np.exp(-806 / ((600 - 273) * temperature + 273))
return arrh_expr
#Calculate time derivative at specific spatial points
def rhsC(c0, c1, c2, T1, dz, p):
return (p[0] * (c2 - 2 * c1 + c0) / dz**2 - p[1] * (c2 - c0) / (2 * dz) - p[2] * k(T1) * c1**2)
def rhsT(T0, T1, T2, c1, dz, p):
return (0.01 * (T2 - 2 * T1 + T0) / dz**2 - 0.5 * (T2 - T0) / (2 * dz) + 2.45 * k(T1) * c1**2)
#Time-stepping loop
for it in range(nt):
c_results[it, :] = C
T_results[it, :] = T
for iz in range(1, nz - 1):
# Calculate concentration and temperature derivatives
dc_dt = rhsC(C[iz - 1], C[iz], C[iz + 1], T[iz], dz, p)
dT_dt = rhsT(T[iz - 1], T[iz], T[iz + 1], C[iz], dz, p)
#Update concentration and temperature using forward Euler
C[iz] += dz * dc_dt
T[iz] += dz * dT_dt
Below is a snippet of my attempt in code B :
class PFR_equation2D(PDE):
name = "PFR_equation2D"
def __init__(self):
# z coordinate
x = Symbol("x")
#time
t = Symbol("t")
input_variables = {"x": x, "t": t}
c = Function("c")(*input_variables)
Temp = Function("Temp")(*input_variables)
def k(Temp):
arrh_exp = 4.68 *exp(-806 / ((600 - 273) * Temp + 273))
return arrh_exp
self.equations = {}
self.equations["PFR_equation_c"] = c.diff(t,1) - 0.01*(c.diff(x,2)) + 0.5*(c.diff(x)) + 1.0*(c**2)
self.equations["PFR_equation_T"] = Temp.diff(t,1) - 0.01*(Temp.diff(x,2)) + 0.5* (Temp.diff(x)) - 2.45*k(Temp)*(c**2)