Value error problem

[AudreyDemafo]

Hello dear community, I am trying to solve a coupled PDE system with deepxde.
My code is the following

import numpy as np
import deepxde as dde
from matplotlib import pyplot as plt
from matplotlib import cm

rho_a = 1.293 # air density
mu_a = 18.6e-6 # Dynamic viscosity of air
cp_a = 1005 # (J/KgK)
rho_v = 0.013 # (kgm-3) 
cp_v = 1996 # (J/KgK)
ka = 0.0285 # air conductivity
beta = 0.973 # Vaporization coefficient of water
R = 0.5 # Sample radius
Deff_a = 0.034 #  Effective diffusivity of air
h1 = 5
a = 1
d = 1
Re = 1

# Define the partial differential model
def NS_equation(x, S):
  u = S[:,0:1]
  v = S[:,1:2]
  w = S[:,2:3]
# , T
# , S[:,3:4]

  # Continuity equation
  gradu_x = dde.grad.jacobian(S, x, i=0, j=0)
  gradv_y = dde.grad.jacobian(S, x, i=1, j=1)
  gradw_z = dde.grad.jacobian(S, x, i=2, j=2)
  continuity = gradu_x + gradv_y + gradw_z

  # Momentum equation
  gradu_y = dde.grad.jacobian(S, x, i=0, j=1)
  gradu_z = dde.grad.jacobian(S, x, i=0, j=2)
  dudt = dde.grad.jacobian(S, x, i=0, j=3)

  gradv_x = dde.grad.jacobian(S, x, i=1, j=0)
  gradv_z = dde.grad.jacobian(S, x, i=1, j=2)
  dvdt = dde.grad.jacobian(S, x, i=1, j=3)

  gradw_x = dde.grad.jacobian(S, x, i=2, j=0)
  gradw_y = dde.grad.jacobian(S, x, i=2, j=1)
  dwdt = dde.grad.jacobian(S, x, i=2, j=3)

  u_hessian_x = dde.grad.hessian(S, x, component=0, i=0, j=0)
  u_hessian_y = dde.grad.hessian(S, x, component=0, i=0, j=1)
  u_hessian_z = dde.grad.hessian(S, x, component=0, i=0, j=2)

  v_hessian_x = dde.grad.hessian(S, x, component=1, i=1, j=0)
  v_hessian_y = dde.grad.hessian(S, x, component=1, i=1, j=1)
  v_hessian_z = dde.grad.hessian(S, x, component=1, i=1, j=2)

  w_hessian_x = dde.grad.hessian(S, x, component=2, i=2, j=0) 
  w_hessian_y = dde.grad.hessian(S, x, component=2, i=2, j=1)
  w_hessian_z = dde.grad.hessian(S, x, component=2, i=2, j=2)

  momentum_x = rho_a*(dudt + u*gradu_x + v*gradu_y + w*gradu_z) - mu_a*(u_hessian_x + u_hessian_y + u_hessian_z)
  momentum_y = rho_a*(dvdt + u*gradv_x + v*gradv_y + w*gradv_z) - mu_a*(v_hessian_x + v_hessian_y + v_hessian_z)
  momentum_z = rho_a*(dwdt + u*gradw_x + v*gradw_y + w*gradw_z) - mu_a*(w_hessian_x + w_hessian_y + w_hessian_z)

  # Energy equation
  gradT_x = dde.grad.jacobian(S, x, i=3, j=0)
  gradT_y = dde.grad.jacobian(S, x, i=3, j=1)
  gradT_z = dde.grad.jacobian(S, x, i=3, j=2)
  dTdt = dde.grad.jacobian(S, x, i=3, j=3)

  grad_gradT_x = ka*dde.grad.hessian(S, x, component=3, i=3, j=0)     
  grad_gradT_y = ka*dde.grad.hessian(S, x, component=3, i=3, j=1)  
  grad_gradT_z = ka*dde.grad.hessian(S, x, component=3, i=3, j=2)   
  diffusion_term = grad_gradT_x + grad_gradT_y + grad_gradT_z
  advection_term = u*gradT_x + v*gradT_y + w*gradT_z
  heat_equation = (rho_a*cp_a + rho_v*cp_v)*dTdt + (rho_a*cp_a + rho_v*cp_v)*advection_term - diffusion_term 

  # # Air humidity dynamics
  gradXa_x = dde.grad.jacobian(S, x, i=4, j=0)
  gradXa_y = dde.grad.jacobian(S, x, i=4, j=1)
  gradXa_z = dde.grad.jacobian(S, x, i=4, j=2)
  dXadt = dde.grad.jacobian(S, x, i=4, j=3)

  grad_gradXa_x = Deff_a*dde.grad.hessian(S, x, component=4, i=4, j=0)
  grad_gradXa_y = Deff_a*dde.grad.hessian(S, x, component=4, i=4, j=1)
  grad_gradXa_z = Deff_a*dde.grad.hessian(S, x, component=4, i=4, j=2)
  diffusion_term1 = grad_gradXa_x + grad_gradXa_y + grad_gradXa_z
  advection_term1 = u*gradXa_x + v*gradXa_y + w*gradXa_z
  Air_humidity = rho_a*dXadt + rho_a*advection_term1 - diffusion_term1

  # # Heat and mass transfer within the product
  gradTp_x = dde.grad.jacobian(S, x, i=5, j=0)
  gradTp_y = dde.grad.jacobian(S, x, i=5, j=1)
  gradTp_z = dde.grad.jacobian(S, x, i=5, j=2)
  dTpdt = dde.grad.jacobian(S, x, i=5, j=3)

  kp = 0.028 + 0.1321*S[:,6:7] + 3.4981e-3*S[:,5:6] - 0.128*S[:,6:7]**2 - 2.24e-5*S[:,5:6]**2
  rho_p = 842.05 - 8.4321*S[:, 6:7]
  cp_p = 0.7494 + 1.228*S[:, 6:7] + 0.04857*S[:,5:6] - 0.7344*S[:,6:7]**2 + 1.0057e-4*S[:,5:6]**2
  Dth_p = ( 821.39 - 821.55*S[:, 6:7] + 990.66*S[:, 6:7]**2 )/rho_p 
  porosity_p = 40.69 + 0.1915*S[:,6:7]
  Deff_p = Dth_p/porosity_p

  grad_gradTp_x = kp*dde.grad.hessian(S, x, component=5, i=5, j=0)
  grad_gradTp_y = kp*dde.grad.hessian(S, x, component=5, i=5, j=1)
  grad_gradTp_z = kp*dde.grad.hessian(S, x, component=5, i=5, j=2)
  diffusion_term2 = grad_gradTp_x + grad_gradTp_y + grad_gradTp_z
  heat_equation2 = (rho_p*cp_p)*dTpdt - diffusion_term2

  gradXp_x = dde.grad.jacobian(S, x, i=6, j=0)
  gradXp_y = dde.grad.jacobian(S, x, i=6, j=1)
  gradXp_z = dde.grad.jacobian(S, x, i=6, j=2)
  dXpdt = dde.grad.jacobian(S, x, i=6, j=3)

  grad_gradXp_x = Deff_p*dde.grad.hessian(S, x, component=6, i=6, j=0)
  grad_gradXp_y = Deff_p*dde.grad.hessian(S, x, component=6, i=6, j=1)
  grad_gradXp_z = Deff_p*dde.grad.hessian(S, x, component=6, i=6, j=2)
  diffusion_term3 = grad_gradXp_x + grad_gradXp_y + grad_gradXp_z
  Moisture_content = (rho_p)*dXpdt - diffusion_term3 

  return [continuity, momentum_x, momentum_y, momentum_z, heat_equation, Air_humidity, heat_equation2, Moisture_content]
  
# The domain
geom1 = dde.geometry.geometry_3d.Sphere(center=(0,0,0), radius=1)
geom3 = dde.geometry.geometry_3d.Cuboid(xmin=(-1, -1, -1), xmax=(2, 2, 2))
geom_CSGD = dde.geometry.CSGDifference(geom3, geom1)
TimeDomain = dde.geometry.TimeDomain(0,2)
geomtime = dde.geometry.GeometryXTime(geom_diff, TimeDomain)

 def u_initial(x):
  return 0.5

def v_initial(x):
  return 0.5

def w_initial(x):
  return 0.5

def Tp_initial(x):
  return 315

def Ta_initial(x):
  return 298.15

def Xa_initial(x):
  return 0.89

def Xp_initial(x):
  return 0.97

# Define the boundary

def boundary_product(x, on_boundary):  # Boundary of the product
  return on_boundary and geom1.on_boundary(x[:3])

def initial_condition_cube(x, on_initial):
  return on_initial 

# Apply the equations to the boundaries
def func_Xp(S):
  rho_p = 842.05 - 8.4321*S[:,6:7]
  Dth_p = ( 821.39 - 821.55*S[:,6:7] + 990.66*S[:,6:7]**2 )/rho_p 
  porosity_p = 40.69 + 0.1915*S[:,6:7]
  Deff_p = Dth_p/porosity_p
  return beta/(rho_p*Deff_p)*(S[:,4:5] - S[:,6:7])

def func_Tp(S):
  kp = 0.028 + 0.1321*S[:,6:7] + 3.4981e-3*S[:,5:6] - 0.128*S[:,6:7]**2 - 2.24e-5*S[:,5:6]**2
  return (h1*np.pi*R**2)*(S[:,3:4] - S[:,5:6])/kp

# Boundary conditions
bc1 = dde.RobinBC(geomtime, func_Xp, boundary_product, component=6) # On the product
bc2 = dde.RobinBC(geomtime, func_Tp, boundary_product, component=5) # On the product
ic_Tp= dde.IC(geomtime, Tp_initial, initial_ellipse,component=5)
ic_Xp = dde.IC(geomtime, Xp_initial, initial_ellipse, component=6)
ic_u = dde.IC(geomtime, u_initial, initial_condition_cube, component=0)
ic_v = dde.IC(geomtime, v_initial, initial_condition_cube, component=1)
ic_w = dde.IC(geomtime, w_initial, initial_condition_cube, component=2)
ic_Ta = dde.IC(geomtime, Ta_initial, initial_condition_cube, component=3)
ic_Xa = dde.IC(geomtime, Xa_initial, initial_condition_cube, component=4)

#bc1,bc2,,
data = dde.data.TimePDE(geomtime,
                        NS_equation,
                        [bc1,bc2,ic_Tp,ic_Xp,ic_u,ic_v,ic_w,ic_Ta,ic_Xa],
                        num_domain = 5000,
                        num_boundary = 2000,
                        num_initial=20000,
                        num_test=1000
                    )

nnet = dde.maps.FNN([4] + [100]*6 + [7], "relu", "Glorot uniform")

model = dde.Model(data, nnet)

model.compile("adam", lr = 0.0001, loss="mse") 

model.train(iterations=20, display_every=1)

losshistory, train_state = model.train(iterations=20)

dde.saveplot(
    losshistory, train_state, issave=True, isplot=True)

And I am getting the following error

ValueError: i=4 is not valid.

I have checked my pde function and everything seems correct to me. I don't know how to fix the error. Any help from you will be appreciated.
@lululxvi, @pescap, please kindly suggest. Thank you!

[lululxvi]

i=4 means the 5th component. Your network only has 4 inputs.

[MosDefMan]

I get for my following code the same error:

import deepxde as dde
import numpy as np
import tensorflow as tf

# Define the competitive Langmuir isotherm parameters and solid phase
def pde(x, y):
    C1, q1 = y[:, 0:1], y[:, 1:2]
    C2, q2 = y[:, 2:3], y[:, 3:4]

    C1_t, C1_xx = dde.grad.jacobian(y, x, i=0, j=1), dde.grad.hessian(y, x, i=0, j=0)
    q1_t = dde.grad.jacobian(y, x, i=1, j=1)
    C2_t, C2_xx = dde.grad.jacobian(y, x, i=2, j=1), dde.grad.hessian(y, x, i=2, j=0)
    q2_t = dde.grad.jacobian(y, x, i=3, j=1)

    denom = 1 + K_L1 * C1 + K_L2 * C2

    q1_eq = q_m1 * K_L1 * C1 / denom
    q2_eq = q_m2 * K_L2 * C2 / denom

    return [
        C1_t - D1 * C1_xx + k_ad * (q1 - q1_eq),
        q1_t - k_ad * (q1_eq - q1),
        C2_t - D2 * C2_xx + k_ad * (q2 - q2_eq),
        q2_t - k_ad * (q2_eq - q2),
    ]

# Initial conditions
def ic_C1(x):
    return np.sin(np.pi * x[:, 0:1])
def ic_C2(x):
    return np.sin(2 * np.pi * x[:, 0:1])
def ic_q1(x):
    return np.zeros_like(x[:, 0:1])
def ic_q2(x):
    return np.zeros_like(x[:, 0:1])

# Constant boundary conditions
def bc_left_C1(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)

def bc_right_C1(x, on_boundary):
    return on_boundary and np.isclose(x[0], 1)

def bc_left_C2(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)

def bc_right_C2(x, on_boundary):
    return on_boundary and np.isclose(x[0], 1)

# Boundary values
C1_left = 1.0
C1_right = 0.0
C2_left = 1.0
C2_right = 0.0

# Define the geometry and time domain
geom = dde.geometry.Interval(0, 1)
timedomain = dde.geometry.TimeDomain(0, 1)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

# Boundary and initial conditions
bc_C1_left = dde.icbc.DirichletBC(geomtime, lambda x: C1_left, bc_left_C1, component=0)
bc_C1_right = dde.icbc.DirichletBC(geomtime, lambda x: C1_right, bc_right_C1, component=0)
bc_C2_left = dde.icbc.DirichletBC(geomtime, lambda x: C2_left, bc_left_C2, component=2)
bc_C2_right = dde.icbc.DirichletBC(geomtime, lambda x: C2_right, bc_right_C2, component=2)

ic_C1 = dde.icbc.IC(geomtime, ic_C1, lambda _, on_initial: on_initial, component=0)
ic_q1 = dde.icbc.IC(geomtime, ic_q1, lambda _, on_initial: on_initial, component=1)
ic_C2 = dde.icbc.IC(geomtime, ic_C2, lambda _, on_initial: on_initial, component=2)
ic_q2 = dde.icbc.IC(geomtime, ic_q2, lambda _, on_initial: on_initial, component=3)

# Measured data points (example data)
data_points = np.array([[0.2, 0.1], [0.5, 0.2], [0.7, 0.3]])
data_values = np.random.rand(data_points.shape[0], 4)

observe_C1 = dde.icbc.PointSetBC(data_points, data_values[:, 0:1], component=0)
observe_q1 = dde.icbc.PointSetBC(data_points, data_values[:, 1:2], component=1)
observe_C2 = dde.icbc.PointSetBC(data_points, data_values[:, 2:3], component=2)
observe_q2 = dde.icbc.PointSetBC(data_points, data_values[:, 3:4], component=3)

# Combine the boundary conditions
data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc_C1_left, bc_C1_right, bc_C2_left, bc_C2_right, ic_C1, ic_q1, ic_C2, ic_q2, observe_C1, observe_q1, observe_C2, observe_q2],
    num_domain=400,
    num_boundary=2,
    num_initial=50,
)

# Define the neural network
net = dde.maps.FNN([2] + [64] * 3 + [4], "tanh", "Glorot uniform")

# Define the unknown parameters q_m, K_L, and k_ad
q_m1, K_L1 = dde.Variable(1.0), dde.Variable(1.0)
q_m2, K_L2 = dde.Variable(1.0), dde.Variable(1.0)
k_ad = dde.Variable(1.0)

# Diffusion coefficients
D1, D2 = 1.0, 1.0

# Define the model
model = dde.Model(data, net)

# Compile the model
model.compile("adam", lr=1e-3, loss="MSE", external_trainable_variables=[q_m1, K_L1, q_m2, K_L2, k_ad])
losshistory, train_state = model.train(epochs=5000)

# Switch to L-BFGS-B optimizer
model.compile("L-BFGS-B", loss="MSE", external_trainable_variables=[q_m1, K_L1, q_m2, K_L2, k_ad])
losshistory, train_state = model.train()

# Display the identified parameters
print(f"Identified parameters:")
print(f"Component 1: q_m={q_m1.value()}, K_L={K_L1.value()}")
print(f"Component 2: q_m={q_m2.value()}, K_L={K_L2.value()}")
print(f"k_ad={k_ad.value()}")

# Plot the results
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

# Predict and visualize the results
X = geomtime.random_points(1000)
y = model.predict(X)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X[:, 0], X[:, 1], y[:, 0], label='C1')
ax.scatter(X[:, 0], X[:, 1], y[:, 1], label='q1')
ax.scatter(X[:, 0], X[:, 1], y[:, 2], label='C2')
ax.scatter(X[:, 0], X[:, 1], y[:, 3], label='q2')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('Concentration / Adsorption')
plt.legend()
plt.show()

Error produced:

  File ~\anaconda3\envs\sioux\lib\site-packages\deepxde\gradients\jacobian.py:59 in __call__
    raise ValueError("j={} is not valid.".format(j))

ValueError: j=2 is not valid.

What am I doing wrong here?
@lululxvi, @pescap, I would appreciate any insight.