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.

Non-linear dynamics tend to allow positive feed-back leading to explosions in finite time. Which gives those overflow errors. This might be a problem of the code, could however also be the behavior of the exact solution that is better approximated with denser grids. What can you say about the exact solution or the expected behavior of the physical original for the model? — Lutz Lehmann, Jan 28 '21 at 12:06
Well, applying the conservation of mass between inlet and outlet i can calculate the velocity at the outlet that is smaller than the inlet one. My Lax-Friedrichs seems is correct in calculating this and i think it works well. The upwind seems to predict a constant value for the velocity and the pressure in the manifold; the Lax-Wendroff doesn't work at all. Because of that i think it's a matter of the code... — SimoneC618, Jan 28 '21 at 14:39

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

0 Answers0