runtime error - Finding RuntimeWarning: overflow encountered in double_scalars while running numerical schemes for fluid dynamic study

Question

i'm writing a code to solve Hyperbolic differential equations with different numerica methods such as Lax-Friederichs, Lax-Wendroff and Upwind scheme. During the calculation i often obtain this type of error:

RuntimeWarning: overflow encountered in double_scalars

that seems to disappear when i reduce the dimensions of matrix. Here i attach my code:

for i in range (0,nt):
#inlet
rho[0,i] = P_inlet/(R*T_inlet)
u[0,i] = u_inlet
P[0,i] = P_inlet
T[0,i] = T_inlet
Ac[0,0] = A_var_list[0]

Q1[0,i] = rho[0,i]
Q2[0,i] = rho[0,i] * u[0,i]
Q3[0,i] = (1/2)*(rho[0,i])*(u[0,i]**2) + (P[0,i]/(k-1))

F1[0,i] = rho[0,i] * u[0,i]
F2[0,i] = (1/2)*(rho[0,i])*(u[0,i]**2) + P[0,i]
F3[0,i] = u[0,i] * ((1/2)*(rho[0,i])*(u[0,i]**2) + (k*P[0,i]/(k-1)))


#outlet
rho[nx-1,i] = rho_outlet
P[nx-1,i] = P_outlet
u[nx-1,i] = u_outlet
T[nx-1,i] = T_outlet
Q1[nx-1,i] = rho[nx-1,i]
Q2[nx-1,i] = rho[nx-1,i]*u[nx-1,i]
Q3[nx-1,i] = (1/2)*rho[nx-1,i]*u[nx-1,i] + (P[nx-1,i]/(k-1))
F1[nx-1,i] = rho[nx-1,i] * u[nx-1,i]
F2[nx-1,i] = (1/2)*rho[nx-1,i]*(u[nx-1,i]**2) + P[nx-1,i]
F3[nx-1,i] = u[nx-1,i] * ((1/2)*(rho[nx-1,i])*(u[nx-1,i]**2) + (k*P[nx-1,i]/(k-1)))

#manifold
for i in range (1,nx-1):
rho[i,0] = P_inlet/(R*Tw[i])
u[i,0] = u_inlet
P[i,0] = P_inlet
Ac[i,0] = A_var_list[i]

Q1[i,0] = rho[i,0]
Q2[i,0] = rho[i,0] * u[i,0]
Q3[i,0] = (1 / 2) * (rho[i,0]) * (u[i,0] ** 2) + (P[i,0] / (k - 1))

F1[i, 0] = rho[i, 0] * u[i, 0]
F2[i, 0] = (1 / 2) * (rho[i, 0]) * (u[i, 0] ** 2) + P[i, 0]
F3[i, 0] = u[i, 0] * ((1 / 2) * (rho[i, 0]) * (u[i, 0] ** 2) + (k * P[i, 0] / (k - 1)))

S1[i, 0] = -rho[i, 0] * u[i, 0] * (Ac[i, 0] - Ac[i - 1, 0])
S2[i, 0] = -(rho[i, 0] * ((u[i, 0] ** 2) / (Ac[i, 0])) * (Ac[i, 0] - Ac[i - 1, 0])) - (
            (frict_fact * np.pi * rho[i, 0] * d[i] * u[i, 0] ** 2) / (2 * Ac[i, 0]))
S3[i, 0] = - (u[i, 0] * (rho[i, 0] * ((u[i, 0] ** 2) / 2) + (k * P[i, 0] / (k - 1))) * (
            (Ac[i, 0] - Ac[i - 1, 0]) / Ac[i, 0])) + (Lambda * np.pi * d[i] * (Tw[i] - T[i, 0]) / Ac[i, 0])
def Upwind():
for n in range (0,nt-1):
    for i in range (1,nx):

        Q1[i,n+1] = Q1[i-1,n]-((F1[i,n] - F1[i-1,n])/Dx)*Dt + (S1[i,n]-S1[i-1,n])*Dt
        Q2[i, n + 1] = Q2[i-1, n] - ((F2[i, n] - F2[i - 1, n]) / Dx) * Dt + (S2[i, n] - S2[i - 1, n]) * Dt
        Q3[i, n + 1] = Q3[i-1, n] - ((F3[i, n] - F3[i - 1, n]) / Dx) * Dt + (S3[i, n] - S3[i - 1, n]) * Dt

        rho[i, n+1] = Q1[i, n+1]
        u[i, n+1] = Q2[i, n+1] / rho[i, n+1]
        P[i, n+1] = (Q3[i, n+1] - 0.5 * rho[i, n+1] * u[i, n+1] ** 2) * (k - 1)
        T[i, n+1] = P[i, n+1] / (R * rho[i, n+1])

        F1[i,n+1] = Q2[i,n+1]
        F2[i,n+1] = rho[i,n+1]*((u[i,n+1]**2)/2) +P[i,n+1]
        F3[i, n + 1] = u[i, n + 1] * (
                    (rho[i, n + 1] * ((u[i, n + 1] ** 2) / 2)) + (k * P[i , n + 1] / (k - 1)))

        S1[i, n + 1] = -rho[i, n + 1] * u[i, n + 1] * (Ac[i, 0] - Ac[i-1, 0])
        S2[i, n + 1] = - (rho[i, n + 1] * (
                (u[i, n + 1] ** 2) / (Ac[i, 0])) * (Ac[i, 0] - Ac[i-1, 0])) - ((
                    (frict_fact * np.pi * rho[i, n + 1] * d[i] * (u[i, n + 1] ** 2)) / (2 * Ac[i, 0])))
        S3[i, n + 1] = -(u[i, n + 1] * (
                    rho[i, n + 1] * ((u[i, n + 1] ** 2) / 2) + (k * P[i, n + 1] / (k - 1))) * (
                                         (Ac[i , 0] - Ac[i-1, 0]) / Ac[i, 0])) + (
                                          Lambda * np.pi * d[i ] * (Tw[i] - T[i, 0]) / Ac[i, 0])



plt.figure(1)
plt.plot(P[:, nt - 1])
plt.figure(2)
plt.plot(u[:, nt - 1])
def Lax_Friedrichs():
for n in range (1,nt):
    for i in range (1,nx-1):
        F1_m1 = 0.5 * (F1[i, n - 1] + F1[i - 1, n - 1])
        F2_m1 = 0.5 * (F2[i, n - 1] + F2[i - 1, n - 1])
        F3_m1 = 0.5 * (F3[i, n - 1] + F3[i - 1, n - 1])

        S1_m1 = 0.5 * (S1[i, 0] + S1[i - 1, 0])
        S2_m1 = 0.5 * (S2[i, 0] + S2[i - 1, 0])
        S3_m1 = 0.5 * (S3[i, 0] + S3[i - 1, 0])


        F1_p1 = 0.5 * (F1[i + 1, n - 1] + F1[i, n - 1])
        F2_p1 = 0.5 * (F2[i + 1, n - 1] + F2[i, n - 1])
        F3_p1 = 0.5 * (F3[i + 1, n - 1] + F3[i, n - 1])

        S1_p1 = 0.5 * (S1[i + 1, n - 1] + S1[i, n - 1])
        S2_p1 = 0.5 * (S2[i + 1, n - 1] + S2[i, n - 1])
        S3_p1 = 0.5 * (S3[i + 1, n - 1] + S3[i, n - 1])


        Q1[i, n] = 0.5 * (Q1[i - 1, n - 1] + Q1[i + 1, n - 1]) - Dt/Dx * (F1_p1 - F1_m1) + (S1_p1 - S1_m1) * Dt
        Q2[i, n] = 0.5 * (Q2[i - 1, n - 1] + Q2[i + 1, n - 1]) - Dt/Dx * (F2_p1 - F2_m1) + (S2_p1 - S2_m1) * Dt
        Q3[i, n] = 0.5 * (Q3[i - 1, n - 1] + Q3[i + 1, n - 1]) - Dt/Dx * (F3_p1 - F3_m1) + (S3_p1 - S3_m1) * Dt

        rho[i, n] = Q1[i, n]
        u[i, n] = Q2[i, n] / rho[i, n]
        P[i, n] = (Q3[i, n] - 0.5 * rho[i, n] * u[i, n] ** 2) * (k - 1)
        T[i, n] = P[i, n] / (R * rho[i, n])

        F1[i, n] = Q2[i, n]
        F2[i, n] = rho[i, n] * ((u[i, n] ** 2) / 2) + P[i, n]
        F3[i, n] = u[i, n] * (
                (rho[i, n] * ((u[i, n] ** 2) / 2)) + (k * P[i, n] / (k - 1)))

        S1[i, n] = -rho[i, n] * u[i, n] * (Ac[i, 0] - Ac[i - 1, 0])
        S2[i, n] = - (rho[i, n] * (
                (u[i, n] ** 2) / (Ac[i, 0])) * (Ac[i, 0] - Ac[i - 1, 0])) - ((
                (frict_fact * np.pi * rho[i, n] * d[i] * (u[i, n] ** 2)) / (2 * Ac[i, 0])))
        S3[i, n] = -(u[i, n] * (
                rho[i, n] * ((u[i, n] ** 2) / 2) + (k * P[i, n] / (k - 1))) * (
                                 (Ac[i, 0] - Ac[i - 1, 0]) / Ac[i, 0])) + (
                               Lambda * np.pi * d[i] * (Tw[i] - T[i, 0]) / Ac[i, 0])

# Plot
plt.figure(1)
plt.plot(P[:, nt - 1])
plt.figure(2)
plt.plot(u[:, nt - 1])

def Lax_Wendroff():
for n in range (0,nt-1):
    for i in range (1,nx-1):

        Q1_plus_half = (1 / 2) * (Q1[i, n] + Q1[i + 1, n]) - (Dt / (2 * Dx)) * (F1[i + 1, n] - F1[i, n]) + (
                    S1[i + 1, n] - S1[i, n]) * Dt
        Q1_less_half = (1 / 2) * (Q1[i, n] + Q1[i - 1, n]) - (Dt / (2 * Dx)) * (F1[i, n] - F1[i - 1, n]) + (
                    S1[i, n] - S1[i - 1, n]) * Dt
        Q2_plus_half = (1 / 2) * (Q2[i-1, n] + Q2[i + 1, n]) - (Dt / (2 * Dx)) * (F2[i + 1, n] - F2[i, n]) + (
                    S2[i + 1, n] - S2[i, n]) * Dt
        Q2_less_half = (1 / 2) * (Q2[i, n] + Q2[i - 1, n]) - (Dt / (2 * Dx)) * (F2[i, n] - F2[i - 1, n]) + (
                    S2[i, n] - S2[i - 1, n]) * Dt
        Q3_plus_half = (1 / 2) * (Q3[i, n] + Q3[i + 1, n]) - (Dt / (2 * Dx)) * (F3[i + 1, n] - F3[i, n]) + (
                    S3[i + 1, n] - S3[i, n]) * Dt
        Q3_less_half = (1 / 2) * (Q3[i, n] + Q3[i - 1, n]) - (Dt / (2 * Dx)) * (F3[i, n] - F3[i - 1, n]) + (
                    S3[i, n] - S3[i - 1, n]) * Dt

        rho_less_half = Q1_less_half
        u_less_half = Q2_less_half / rho_less_half
        P_less_half = (Q3_less_half - ((1 / 2) * rho_less_half * (u_less_half ** 2) / 2)) * (k - 1)

        F1_less_half = rho_less_half * u_less_half
        F2_less_half = rho_less_half * ((u_less_half ** 2) / 2) + P_less_half
        F3_less_half = u_less_half * ((rho_less_half * ((u_less_half ** 2) / 2)) + (k * P_less_half / (k - 1)))

        rho_plus_half = Q1_plus_half
        u_plus_half = Q2_plus_half / rho_plus_half
        P_plus_half = (Q3_plus_half - ((1 / 2) * rho_plus_half * (u_plus_half ** 2) / 2)) * (k - 1)

        F1_plus_half = rho_plus_half * u_plus_half
        F2_plus_half = rho_plus_half * ((u_plus_half ** 2) / 2) + P_plus_half
        F3_plus_half = u_plus_half * ((rho_plus_half * ((u_plus_half ** 2) / 2)) + (k * P_plus_half / (k - 1)))

        # I termini sorgente da mettere dentro l'equazione finale di Q li calcolo come medie delle variabili nel condotto

        S1_less_half = 0.5 * (S1[i - 1, n] + S1[i, n])
        S2_less_half = 0.5 * (S2[i - 1, n] + S2[i, n])
        S3_less_half = 0.5 * (S3[i - 1, n] + S3[i, n])

        S1_plus_half = 0.5 * (S1[i + 1, n] + S1[i, n])
        S2_plus_half = 0.5 * (S2[i + 1, n] + S2[i, n])
        S3_plus_half = 0.5 * (S3[i + 1, n] + S3[i, n])

        """S1_less_half = Q1_less_half + F1_less_half
        S2_less_half = Q2_less_half + F2_less_half
        S3_less_half = Q3_less_half + F3_less_half

        S1_plus_half = Q1_plus_half + F1_plus_half
        S2_plus_half = Q2_plus_half + F2_plus_half
        S3_plus_half = Q3_plus_half + F3_plus_half"""

        Q1[i , n + 1] = Q1[i, n] - (Dt / Dx) * (F1_plus_half - F1_less_half) - (S1_plus_half - S1_less_half) * Dt
        Q2[i, n + 1] = Q2[i, n] - (Dt / Dx) * (F2_plus_half - F2_less_half) - (S2_plus_half - S2_less_half) * Dt
        Q3[i, n + 1] = Q3[i, n] - (Dt / Dx) * (F3_plus_half - F3_less_half) - (S3_plus_half - S3_less_half) * Dt

        rho[i, n + 1] = Q1[i, n + 1]
        u[i, n + 1] = Q2[i, n + 1] / rho[i, n + 1]
        P[i, n + 1] = (Q3[i, n + 1] - 0.5 * rho[i, n + 1] * (u[i, n + 1] ** 2)) * (k - 1)

        F1[i, n + 1] = rho[i, n + 1] * u[i, n + 1]
        F2[i, n + 1] = rho[i, n + 1] * ((u[i, n + 1] ** 2) / 2) + P[i, n + 1]
        F3[i, n+1] = u[i, n+1] * (
                (rho[i, n+1] * ((u[i, n+1] ** 2) / 2)) + (k * P[i, n+1] / (k - 1)))

        S1[i, n+1] = -rho[i, n+1] * u[i, n+1] * (Ac[i, 0] - Ac[i - 1, 0])
        S2[i, n+1] = - (rho[i, n+1] * (
                (u[i, n+1] ** 2) / (Ac[i, 0])) * (Ac[i, 0] - Ac[i - 1, 0])) - ((
                (frict_fact * np.pi * rho[i, n+1] * d[i] * (u[i, n+1] ** 2)) / (2 * Ac[i, 0])))
        S3[i, n+1] = -(u[i, n+1] * (
                rho[i, n+1] * ((u[i, n+1] ** 2) / 2) + (k * P[i, n+1] / (k - 1))) * (
                             (Ac[i, 0] - Ac[i - 1, 0]) / Ac[i, 0])) + (
                           Lambda * np.pi * d[i] * (Tw[i] - T[i, 0]) / Ac[i, 0])


# Plot
plt.figure(1)
plt.plot(P[:, nt - 1])
plt.figure(2)
plt.plot(u[:, nt - 1])

I'm pretty sure that's a matter of indices but i havent't found the solution yet. Hope you can help me.

question from:https://stackoverflow.com/questions/65936262/finding-runtimewarning-overflow-encountered-in-double-scalars-while-running-num

Answer 1

Waitting for answers