Solution of PDAE

[boyeac]

Hi @lululxvi

Thank you in advance for your attention and for making this tool available.

I'm trying to solve a complex problem described by a system of partial differential equations and algebraic (PDAE). My difficulty is in getting Deepxde to understand the different initial and boundary conditions that the system requires.

Of the problems I found on git, I haven't seen any algebraic equations that need to be solved together with the partial ones.

Do you believe it is possible to resolve these types of issues with the current version of Deepxde? Are there any examples that you have already managed to solve?

At the moment, I have an error like: "File "...\deepxde\boundary_conditions.py", line 29, in filter return X[self.on_boundary(X, self.geom.on_boundary(X))]

IndexError: arrays used as indices must be integer (or boolean) type".

The code is:

from future import absolute_import from future import division from future import print_function import numpy as np

import deepxde as dde

Functions

def main():

Constants

epsb       = 0.38
Dax        = 5.46e-4
Rp         = 2e-3
dp         = 2*Rp 
apM        = 3/Rp
kf         = 1.32e-2
visc       = 1.48e-5
Mw1        = 44e-1 
Mw2        = 16e-3
Mw3        = 28e-3
lbda       = 14.6
Cpmix      = 37.06
Rg         = 8.31415
hf         = 28.6
hw         = 120
dw         = 0.0211
Cvmix      = 28.74
Dp1        = 4.49e-6
Dp2        = 4.49e-6
Dp3        = 4.49e-6
rhoap      = 850
epsp       = 0.433
Bii1       = (apM*kf*Rp**2)/(15*epsp*Dp1)
Bii2       = (apM*kf*Rp**2)/(15*epsp*Dp2)
Bii3       = (apM*kf*Rp**2)/(15*epsp*Dp3)
Dm10       = 1.24e7
Dm20       = 1.30e7
Dm30       = 1.59e7
Edm1       = 11424
Edm2       = 8042
Edm3       = 6850
kinf1_1    = 1.198e-5 
kinf1_2    = 7.451e-5
kinf1_3    = 6.174e-5
kinf2_1    = 8.28e-5
kinf2_2    = 6.558e-6
kinf2_3    = 2.395e-4
MDH1_1     = 25387
MDH1_2     = 17871
MDH1_3     = 15222
MDH2_1     = 26665
MDH2_2     = 29857
MDH2_3     = 17110
qmax1_1    = 6.023 
qmax1_2    = 4.662
qmax1_3    = 4.141
qmax2_1    = 1.248 
qmax2_2    = 1.291
qmax2_3    = 1.394
roap       = 850
rob        = (1-epsb)*roap
Cv1        = 28.919
Cv2        = 27.663
Cv3        = 20.800
Cps        = 1046
row        = 8238
Cpw        = 500
alfaw      = 456
alfawl     = 496
U          = 40
Tinf       = 303
Phigh      = 1.9e5
Tinlet     = 303
Cinlet     = Phigh/Rg/Tinlet
y1feed     = 0.865
y2feed     = 0.135
y3feed     = 0
QfeedSLPM  = 0.5;
Dw         = 0.0211;
bed_area   = 3.1415*(Dw/2)**2
u0_inlet   = QfeedSLPM*Tinlet/273*1e5/Phigh/1000/60/bed_area

# PDAE
def pde(x, u):

    # Dependent variables
    y1    = u[:,   0:1]
    y2    = u[:,   1:2]
    y3    = u[:,   2:3]
    Cg1   = u[:,   3:4]
    Cg2   = u[:,   4:5]
    Cg3   = u[:,   5:6]
    Cgt   = u[:,   6:7]
    u0    = u[:,   7:8]
    P     = u[:,   8:9]
    Tg    = u[:,  9:10] 
    Tp    = u[:, 10:11] 
    Tw    = u[:, 11:12]
    Cmp1  = u[:, 12:13]
    Cmp2  = u[:, 13:14]
    Cmp3  = u[:, 14:15]
    q1    = u[:, 15:16]
    q2    = u[:, 16:17]
    q3    = u[:, 17:18]

    # Derivatives
    # dyidz - Jacobian
    dy1dz       = dde.grad.jacobian(u, x,  i=0,  j=0)
    dy2dz       = dde.grad.jacobian(u, x,  i=1,  j=0)
    dy3dz       = dde.grad.jacobian(u, x,  i=2,  j=0)
    dCg1dz      = dde.grad.jacobian(u, x,  i=3,  j=0)
    dCg2dz      = dde.grad.jacobian(u, x,  i=4,  j=0)
    dCg3dz      = dde.grad.jacobian(u, x,  i=5,  j=0)
    dCgTdz      = dde.grad.jacobian(u, x,  i=6,  j=0)
    du0dz       = dde.grad.jacobian(u, x,  i=7,  j=0)
    dPdz        = dde.grad.jacobian(u, x,  i=8,  j=0)
    dTgdz       = dde.grad.jacobian(u, x,  i=9,  j=0)
    #dTpdz       = dde.grad.jacobian(u, x,  i=10, j=0)
    dCg1dt      = dde.grad.jacobian(u, x,  i=3,  j=1)
    dCg2dt      = dde.grad.jacobian(u, x,  i=4,  j=1)
    dCg3dt      = dde.grad.jacobian(u, x,  i=5,  j=1)
    dCgTdt      = dde.grad.jacobian(u, x,  i=6,  j=1)
    dTgdt       = dde.grad.jacobian(u, x,  i=9,  j=1)
    dTpdt       = dde.grad.jacobian(u, x,  i=10, j=1)
    dTwdt       = dde.grad.jacobian(u, x,  i=11, j=1)
    dCmp1dt     = dde.grad.jacobian(u, x,  i=12, j=1)
    dCmp2dt     = dde.grad.jacobian(u, x,  i=13, j=1)
    dCmp3dt     = dde.grad.jacobian(u, x,  i=14, j=1)
    dq1dt       = dde.grad.jacobian(u, x,  i=15, j=1)
    dq2dt       = dde.grad.jacobian(u, x,  i=16, j=1)
    dq3dt       = dde.grad.jacobian(u, x,  i=17, j=1)

    # d2yidz2 - Hessian
    d2y1dzz     = dde.grad.hessian(u, x, component=0, i=0, j=0)
    d2y2dzz     = dde.grad.hessian(u, x, component=1, i=0, j=0)
    d2y3dzz     = dde.grad.hessian(u, x, component=2, i=0, j=0)
    d2Tgdzz     = dde.grad.hessian(u, x, component=9, i=0, j=0)

    # Algebraic Equations
    # Total density
    rho = Cgt * (Mw1*y1+Mw2*y2+Mw3*y3)   
    # 
    Cs1 = Bii1*Cg1+Cmp1/(1+Bii1)
    Cs2 = Bii2*Cg2+Cmp2/(1+Bii2)
    Cs3 = Bii3*Cg3+Cmp3/(1+Bii3)
    # 
    Dm1 = Dm10*np.exp(-Edm1/Rg/Tp)
    Dm2 = Dm20*np.exp(-Edm2/Rg/Tp)
    Dm3 = Dm30*np.exp(-Edm3/Rg/Tp)

    # 
    Pmpbar1 = Cmp1*Rg*Tp/1e5
    Pmpbar2 = Cmp2*Rg*Tp/1e5
    Pmpbar3 = Cmp3*Rg*Tp/1e5
    # 
    keq1_1 = kinf1_1*np.exp(MDH1_1/Rg/Tp)
    keq1_2 = kinf1_2*np.exp(MDH1_2/Rg/Tp)
    keq1_3 = kinf1_3*np.exp(MDH1_3/Rg/Tp)
    # 
    keq2_1 = kinf2_1*np.exp(MDH2_1/Rg/Tp)
    keq2_2 = kinf2_2*np.exp(MDH2_2/Rg/Tp)
    keq2_3 = kinf2_3*np.exp(MDH2_3/Rg/Tp)
    # 
    qast1 = (qmax1_1*keq1_1*Pmpbar1/(1+keq1_1*Pmpbar1+keq1_2*Pmpbar2+keq1_3*Pmpbar3) 
             + qmax2_1*keq2_1*Pmpbar1/(1+keq2_1*Pmpbar1+keq2_2*Pmpbar2+keq2_3*Pmpbar3)) 

    qast2 = (qmax1_2*keq1_2*Pmpbar2/(1+keq1_1*Pmpbar1+keq1_2*Pmpbar2+keq1_3*Pmpbar3) 
             + qmax2_2*keq2_2*Pmpbar1/(1+keq2_1*Pmpbar1+keq2_2*Pmpbar2+keq2_3*Pmpbar3))

    qast3 = (qmax1_3*keq1_3*Pmpbar3/(1+keq1_1*Pmpbar1+keq1_2*Pmpbar2+keq1_3*Pmpbar3) 
             + qmax2_3*keq2_3*Pmpbar3/(1+keq2_1*Pmpbar1+keq2_2*Pmpbar2+keq2_3*Pmpbar3))

    # Ideal gas law
    P = Cgt*Tg*Rg

    ## Differential equations
    # Component 1 balance
    BC_1 = (epsb*Dax*Cgt*d2y1dzz + dy1dz*dCgTdz 
                            - (u0*dCg1dz + du0dz*Cg1) 
                            - epsb*dCg1dt - (1-epsb)*apM*kf*(Cg1-Cs1))
    # Component 2 balance
    BC_2 =  (epsb*Dax*Cgt*d2y2dzz + dy2dz*dCgTdz 
                             - (u0*dCg2dz + du0dz*Cg2) 
                             - epsb*dCg2dt - (1-epsb)*apM*kf*(Cg2-Cs2))
    # Component 3 balance
    BC_3 = (epsb*Dax*Cgt*d2y3dzz + dy3dz*dCgTdz 
                             - (u0*dCg3dz + du0dz*Cg3) 
                             - epsb*dCg3dt - (1-epsb)*apM*kf*(Cg3-Cs3))
    # Ergun Equation
    Ergun = (dPdz + 150*visc*((1-epsb)**2)*u0/((epsb**3)*(dp**2))
             + 1.75*(1-epsb)*(rho)*abs(u0)*u0/((epsb**3)*dp))
    # Energy balance in the gas phase
    BEFG = (lbda*d2Tgdzz - u0*Cgt*Cpmix*dTgdz
                                + epsb*Rg*Tg*dCgTdt 
                                - (1-epsb)*apM*hf*(Tg-Tp)
                                - (4*hw/dw)*(Tg-Tw)
                                - epsb*Cgt*Cvmix*dTgdt)
    # Component 1 balance in solid phase
    BMFS_1 = (dCmp1dt - (15*Dp1/(Rp**2))*(Cs1-Cmp1)
                                   -(rhoap/epsp)*dq1dt)
    # Component 2 balance in solid phase
    BMFS_2 = (dCmp2dt - (15*Dp2/(Rp**2))*(Cs2-Cmp2)
                                   -(rhoap/epsp)*dq2dt)
    # Component 3 balance in solid phase
    BMFS_3 = (dCmp3dt - (15*Dp3/(Rp**2))*(Cs3-Cmp3)
                                   -(rhoap/epsp)*dq3dt)
    # Material solid phase micro Baalanço
    BMC_1 = (dq1dt - 15*Dm1*(qast1-q1))
    # Material solid phase micro Baalanço
    BMC_2 = (dq2dt - 15*Dm2*(qast2-q2))
    # Material solid phase micro Baalanço
    BMC_3 = (dq3dt - 15*Dm3*(qast3-q3))

    # Energy equation
    EEnergy = ( (1-epsb)*(epsp*(Cmp1*Cv1+Cmp2*Cv2+Cmp3*Cv3) + roap*(q1*Cv1+q2*Cv2+q3*Cv3)+  roap*Cps)*dTpdt 
                       - (1-epsb)*epsp*Rg*Tp*(dCmp1dt+dCmp2dt+dCmp3dt) 
                       - rob*(MDH1_1*dq1dt+MDH1_2*dq2dt+MDH1_3*dq3dt)
                       - (1-epsb)*apM*hf*(Tg-Tp) )
    # Energy balance on the wall
    BEWall = (row*Cpw*dTwdt - alfaw*hw*(Tg-Tw) + alfawl*U*(Tw-Tinf) 
                                 )

    return [BC_1, BC_2, BC_3, 
            Ergun, BEFG,BMFS_1,
            BMFS_2,BMFS_3,
            BMC_1,BMC_2,
            BMC_3,EEnergy, BEWall]

#def Cond_inicialq2(x,u):
#    return dq2dt

def cc_Cgt(x):
    return (x[:, 3:4] + x[:, 4:5] + x[:, 5:6])

def cc_dCg1dz_0(x):

    return (((y1feed * Cinlet) - (x[:, 3:4])) / (-(epsb * Dax)/ u0_inlet))

def cc_dCg2dz_0(x):
    return ((y2feed * Cinlet - (x[:, 4:5])) / (-(epsb * Dax)/ u0_inlet))

def cc_dCg3dz_0(x):
    return ((y3feed * Cinlet - (x[:, 5:6])) / (-(epsb * Dax)/ u0_inlet))

def cc_y1(x):
    return (x[:, 3:4]/x[:, 6:7])

def cc_y2(x):
    return (x[:, 4:5]/x[:, 6:7])

def cc_y3(x):
    return (x[:, 5:6]/x[:, 6:7])

def cc_u0_0(x):
    return ((u0_inlet*Cinlet)/ x[:, 6:7])

def cc_dTgdz_0(x):
    return ((Tinlet - x[:, 9:10]) / (-lbda/(u0_inlet*Cinlet*Cpmix)))

def cc_P(x):
    return (x[:, 6:7]*x[:, 9:10]*Rg)

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

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

def boundary_initial(x, _):
  return  np.isclose(x[0], 0 )

def boundary(x, on_boundary):
  return on_boundary

# Definition of the geometry domain
spatial_domain = dde.geometry.Interval(0, 0.5)
# Time domain
temporal_domain = dde.geometry.TimeDomain(0, 100)
# Spatio temporal domain
spatio_temporal_domain = dde.geometry.GeometryXTime(spatial_domain, temporal_domain)

# Boundary conditions

c_Cgt_0    = dde.DirichletBC(spatio_temporal_domain, cc_Cgt, lambda _, on_boundary: on_boundary_0, component=6)          #Ok
c_Cgt_L    = dde.DirichletBC(spatio_temporal_domain, cc_Cgt, lambda _, on_boundary: on_boundary_L, component=6)          #Ok
c_Cg1_0    = dde.RobinBC(spatio_temporal_domain, cc_dCg1dz_0, lambda _, on_boundary: on_boundary_0, component=3)         #Ok
c_Cg2_0    = dde.RobinBC(spatio_temporal_domain, cc_dCg2dz_0, lambda _, on_boundary: on_boundary_0, component=4)         #Ok
c_Cg3_0    = dde.RobinBC(spatio_temporal_domain, cc_dCg3dz_0, lambda _, on_boundary: on_boundary_0, component=5)         #Ok
c_Cg1_L    = dde.NeumannBC(spatio_temporal_domain, lambda x: 0, lambda _, on_boundary: on_boundary_L, component=3)       #Ok
c_Cg2_L    = dde.NeumannBC(spatio_temporal_domain, lambda x: 0, lambda _, on_boundary: on_boundary_L, component=4)       #Ok
c_Cg3_L    = dde.NeumannBC(spatio_temporal_domain, lambda x: 0, lambda _, on_boundary: on_boundary_L, component=5)       #Ok
c_Tg_0     = dde.RobinBC(spatio_temporal_domain, cc_dTgdz_0, lambda _, on_boundary: on_boundary_0,  component=9)         #Ok
c_Tg_L     = dde.NeumannBC(spatio_temporal_domain, lambda x: 0, lambda _, on_boundary: on_boundary_L, component=9)       #Ok
c_u0_0     = dde.DirichletBC(spatio_temporal_domain, cc_u0_0, lambda _, on_boundary: on_boundary_0, component=7)         #Ok 
c_u0_L     = dde.DirichletBC(spatio_temporal_domain, lambda x: 0, lambda _, on_boundary: on_boundary_L, component=7)    #Ok
c_P_0      = dde.DirichletBC(spatio_temporal_domain, cc_P, lambda _, on_boundary: on_boundary_0, component=8)          #Ok
c_P_L      = dde.DirichletBC(spatio_temporal_domain, cc_P, lambda _, on_boundary: on_boundary_L, component=8) #Ok
c_y1_0     = dde.DirichletBC(spatio_temporal_domain, cc_y1, lambda _, on_boundary: on_boundary_0, component=0)           #Ok
c_y2_0     = dde.DirichletBC(spatio_temporal_domain, cc_y2, lambda _, on_boundary: on_boundary_0, component=1)           #Ok
c_y3_0     = dde.DirichletBC(spatio_temporal_domain, cc_y3, lambda _, on_boundary: on_boundary_0, component=2)           #Ok
c_y1_L     = dde.DirichletBC(spatio_temporal_domain, cc_y1, lambda _, on_boundary: on_boundary_L, component=0)           #Ok
c_y2_L     = dde.DirichletBC(spatio_temporal_domain, cc_y2, lambda _, on_boundary: on_boundary_L, component=1)           #Ok
c_y3_L     = dde.DirichletBC(spatio_temporal_domain, cc_y3, lambda _, on_boundary: on_boundary_L, component=2)           #Ok

# Initial conditions

ic_y1   = dde.IC(spatio_temporal_domain, lambda x: 0, lambda _, on_initial: on_initial, component=0)               #Ok
ic_y2   = dde.IC(spatio_temporal_domain, lambda x: 1, lambda _, on_initial: on_initial, component=1)               #Ok
ic_y3   = dde.IC(spatio_temporal_domain, lambda x: 0, lambda _, on_initial: on_initial, component=2)               #Ok
ic_Cgt  = dde.IC(spatio_temporal_domain, lambda x: Cinlet, lambda _, on_initial: on_initial, component=6)          #Ok
ic_Tg   = dde.IC(spatio_temporal_domain, lambda x: Tinlet, lambda _, on_initial: on_initial, component=9)          #Ok
ic_Tp   = dde.IC(spatio_temporal_domain, lambda x: Tinlet, lambda _, on_initial: on_initial, component=10)         #Ok
ic_Tw   = dde.IC(spatio_temporal_domain, lambda x: Tinlet, lambda _, on_initial: on_initial, component=11)         #Ok
ic_Cmp1 = dde.IC(spatio_temporal_domain, lambda x: 0, lambda _, on_initial: on_initial, component=12)              #Ok
ic_Cmp2 = dde.IC(spatio_temporal_domain, lambda x: Cinlet, lambda _, on_initial: on_initial, component=13)         #Ok
ic_Cmp3 = dde.IC(spatio_temporal_domain, lambda x: 0, lambda _, on_initial: on_initial, component=14)              #Ok
ic_q1   = dde.IC(spatio_temporal_domain, lambda x: 0, lambda _, on_initial: on_initial, component=15)              #Ok
ic_q2   = dde.NeumannBC(spatio_temporal_domain, lambda x: 0, lambda _, on_initial: boundary_initial, component=16)
ic_q3   = dde.IC(spatio_temporal_domain, lambda x: 0, lambda _, on_initial: on_initial, component=17)              #Ok

# Hyperparameters of the net.

data = dde.data.TimePDE(
    spatio_temporal_domain,
    pde,
    [
        c_Cgt_0,
        c_Cgt_L,
        c_Cg1_0,
        c_Cg2_0,
        c_Cg3_0,
        c_Cg1_L,
        c_Cg2_L,
        c_Cg3_L,
        c_Tg_0,
        c_Tg_L,
        c_u0_0,
        c_u0_L,
        c_P_0,
        c_P_L,
        c_y1_0,
        c_y2_0,
        c_y3_0,
        c_y1_L,
        c_y2_L,
        c_y3_L,
        ic_y1,
        ic_y2,
        ic_y3,
        ic_Cgt,
        ic_Tg,
        ic_Tp,
        ic_Tw,
        ic_Cmp1,
        ic_Cmp2,
        ic_Cmp3,
        ic_q1,
        ic_q2,
        ic_q3,
     ],
     num_domain=100,
     num_boundary=100,
     num_initial=100,
     num_test=100,
)
# Preliminary set
net = dde.maps.FNN([1] + 1 * [5] + [1], "tanh", "Glorot normal")
# Build the model
model = dde.Model(data, net)
# Compile the model
model.compile("adam", lr=1e-3)
# 
model.train(epochs=1)
# 
model.compile("L-BFGS-B")
# 
losshistory, train_state = model.train()

if name == "main": main()

End of the code

Thanks for the help!

[boyeac]

Hey @lululxvi

I found that the solution to the problem is possible. Checking the other problems, I realized that I would need to normalize the variables and make a series of modifications that I've already incorporated because they were really wrong.

Now the code runs without errors referring to DeepXDE, but I'm getting some NaN in the functions in response. I believe they come from some normalization error, parameters, or errors in some equation.

If you can give me any ideas, I would appreciate it.

The output values are (typically): "" INFO:tensorflow:Optimization terminated with: Message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT Objective function value: nan Number of iterations: 2864 Number of functions evaluations: 15001 15002 [8.30e+06, 1.48e+06, 3.42e+02, 2.03e+03, 3.96e+05, 1.79e+06, 3.85e+04, 2.07e+02, 2.94e+05, 3.59e+05, 2.95e+05, 3.49e+05, 1.64e+05, 3.30e+04, 1.87e+05, 3.17e+00, nan, nan, 1.94e+07, 5.04e+05, 6.75e+03, 2.78e-06, 1.83e-05, 7.11e-06, nan, nan, nan, 1.03e+05, 9.67e-07, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan] [8.21e+06, 1.47e+06, 3.39e+02, 2.00e+03, 2.35e+07, 2.40e+06, 3.93e+04, 1.56e+03, 4.97e+12, 7.89e+11, 1.66e+08, 1.99e+13, 5.53e+10, 3.27e+04, 6.26e+07, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00] []

Best model at step 0: train loss: inf test loss: inf test metric: 'train' took 425.262300 s ""

Thank you very much.

[lululxvi]

Your loss like 8.30e+06 is too large, so the SGD may just blow up. You can either scale your problem such that the loss in step 1 is O(1), or you can use a small loss_weights in Model.compile().

[zyz-rise]

@lululxvi How can I scale the components in the PDE to the same order of magnitude? For example, O(1)

[lululxvi]

See FAQ: Q: I failed to train the network or get the right solution, e.g., large training loss, unbalanced losses.