Euler eqs. example

[lkampoli]

Hello,
I'd like to ask whether would be possible to add a simple example for Euler equations, even 1D, for example a nozzle flow or 1d shock flow, a riemann problem, please?
Or some hints on how to do it?

Kind regards,
Lorenzo Campoli

[lululxvi]

You can check the examples of Burgers, Poisson, diffusion equation, etc. The Euler equation is very similar to them. Let me know if you have any difficulty.

[lkampoli]

Hello, thanks for your fast response and help.
I'm trying to setup the 1D Euler shock wave case as presented in https://doi.org/10.1016/j.cma.2019.112789

At the moment, I have this, which doesn't really work. Am I on the good way?

from future import absolute_import
from future import division
from future import print_function

import numpy as np

import sys
sys.path.insert(0, '..')

import deepxde as dde
from deepxde.backend import tf

def main():
def Euler_system(x, y):
"""1D Euler equations system.
d(rho)/dt + d(rhou)/dt = 0
d(rhou)/dt + d(rhouu + p)/dx = 0
d(rhoE)/dt + d(u(rhoe + p)/dx = 0
p = (gamma - 1) * (rho * E - 0.5 * |u|^2)
"""
P = y[:,0:1]
rho = y[:,1:2]
u = y[:,2:3]
E = y[:,3:4]

    rho_x = tf.gradients(rho, x)[0]
    drho_x, drho_t = rho_x[:, 0:1], rho_x[:, 1:2]

    rhou_x = tf.gradients(rho*u, x)[0]
    drhou_x, drhou_t = rhou[:, 0:1], rhou[:, 1:2]

    rhoE_x = tf.gradients(rho*E, x)[0]
    drhoE_x, drhoE_t = rhoE[:, 0:1], rhoE[:, 1:2]

    mass_flow_grad  = tf.gradients(rho*u, x)[0]
    momentum_grad   = tf.gradients((rho*u*u + P), x)[0]
    energy_grad     = tf.gradients((rho*E + P)*u, x)[0]

    gamma = 1.4
    state_res = P - rho*(gamma-1)*(E-0.5*gamma*u*u)

    return [drho_t  + mass_flow_grad, 
            drhou_t + momentum_grad,
            drhoE_t + energy_grad,
            state_res
    ]

def boundary(_, on_initial):
    return on_initial

def func(x):
    """Initial solution.
    rho = 1.4 if x < 0.5
    rho = 1.0 if x > 0.5
    u   = 0.1
    P   = 1.0
    E   = 1/rho * P/(gamma-1) + 0.5 * u * u
    """
    P = 1.0
    if x<0.5:
        rho = 1.4
    else:
        rho = 1.0
    u = 0.1
    gamma = 1.4
    E = 1/rho * P/(gamma-1) + 0.5 * u * u
    return [P*np.ones(len(x),1),
            rho*np.ones(len(x),1),
            u*np.ones(len(x),1),
            E*np.ones(len(x),1)
    ]

def solution(x):
    """Exact solution.
    rho = 1.4 if x < 0.5 + 0.1 *t
    rho = 1.0 if x > 0.5 + 0.1 *t
    u   = 0.1
    P   = 1.0
    E   = 1/rho * P/(gamma-1) + 0.5 * u * u
    """
    x, t = x[:, 0:1], x[:, 1:]
    P = 1.0
    if x<0.5+0.1*t:
        rho = 1.4
    else:
        rho = 1.0
    u = 0.1
    gamma = 1.4
    E = 1/rho * P/(gamma-1) + 0.5 * u * u
    #return [P, rho, u, E]
    return [P*np.ones(len(x),1),
            rho*np.ones(len(x),1),
            u*np.ones(len(x),1),
            E*np.ones(len(x),1)
    ]

geom = dde.geometry.Interval(0, 1)
timedomain = dde.geometry.TimeDomain(0, 2)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

ic1 = dde.IC(geomtime, func, lambda _, on_initial: on_initial, component=0)
ic2 = dde.IC(geomtime, func, lambda _, on_initial: on_initial, component=1)
ic3 = dde.IC(geomtime, func, lambda _, on_initial: on_initial, component=2)
ic4 = dde.IC(geomtime, func, lambda _, on_initial: on_initial, component=3)

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

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

bc_l1 = dde.DirichletBC(geomtime, solution, lambda _, on_boundary: on_boundary, component=0)
bc_l2 = dde.DirichletBC(geomtime, solution, lambda _, on_boundary: on_boundary, component=1)
bc_l3 = dde.DirichletBC(geomtime, solution, lambda _, on_boundary: on_boundary, component=2)
bc_l4 = dde.DirichletBC(geomtime, solution, lambda _, on_boundary: on_boundary, component=3)
bc_r1 = dde.DirichletBC(geomtime, solution, lambda _, on_boundary: on_boundary, component=0)
bc_r2 = dde.DirichletBC(geomtime, solution, lambda _, on_boundary: on_boundary, component=1)
bc_r3 = dde.DirichletBC(geomtime, solution, lambda _, on_boundary: on_boundary, component=2)
bc_r4 = dde.DirichletBC(geomtime, solution, lambda _, on_boundary: on_boundary, component=3)

data = dde.data.TimePDE(
    geomtime,
    Euler_system,
    [bc_l1, bc_l2, bc_l3, bc_l4, bc_r1, bc_r2, bc_r3, bc_r4, ic1, ic2, ic3, ic4],
    num_domain=400,
    num_boundary=40,
    num_initial=100,
    num_test=10000,
)

layer_size = [1] + [20] * 7 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)
model.compile("adam", lr=0.001, metrics=["l2 relative error"])
losshistory, train_state = model.train(epochs=10000)

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

checkpointer = dde.callbacks.ModelCheckpoint(
    "./model/model.ckpt", verbose=1, save_better_only=True
)
movie = dde.callbacks.MovieDumper(
    "model/movie", [-1], [1], period=100, save_spectrum=True, y_reference=func
)

# Plot PDE residue
x = geom.uniform_points(1000, True)
y = model.predict(x, operator=pde)
plt.figure()
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("PDE residue")
plt.show()

if name == "main":
main()

[lululxvi]

Some suggestions:

• rho, u, p, E are not independent. The network should have 3 outputs, e.g., rho, u, p. Then E can be calculated from them.
• Try a smooth solution first, e.g., the example in Appendix A.
• You need to define a separate func for each IC, see the example Lorenz_inverse.py
• Similarly, each BC should have a separate function. Check the examples with multiple BCs, e.g., Poisson_Neumann_1d.py.

[lkampoli]

Hello, thanks for suggestions. If I understand correctly what you mean, I modified the code as follows, but I get an error about the BC:

File "../deepxde/boundary_conditions.py", line 54, in error
"DirichletBC should output 1D values. Use argument 'component' for different components."
RuntimeError: DirichletBC should output 1D values. Use argument 'component' for different components.

Thank you,
Lorenzo

from future import absolute_import
from future import division
from future import print_function

import numpy as np

import sys
sys.path.insert(0, '..')

import deepxde as dde
from deepxde.backend import tf

def main():
def Euler_system(x, y):
"""1D Euler equations system.
d(rho)/dt + d(rhou)/dt = 0
d(rhou)/dt + d(rhouu + p)/dx = 0
d(rhoE)/dt + d(u(rhoe + p)/dx = 0
p = (gamma - 1) * (rho * E - 0.5 * |u|^2)
"""
rho = y[:, 0:1]
u = y[:, 1:2]
p = y[:, 2:3]

    gamma = 1.4
    E = 1/rho * p/(gamma-1) + 0.5 * u * u

    # First order derivatives
    rho_x = tf.gradients(rho, x)[0]
    drho_x, drho_t = rho_x[:, 0:1], rho_x[:, 1:2]

    rhou_x = tf.gradients(rho*u, x)[0]
    drhou_x, drhou_t = rhou_x[:, 0:1], rhou_x[:, 1:2]

    rhoE_x = tf.gradients(rho*E, x)[0]
    drhoE_x, drhoE_t = rhoE_x[:, 0:1], rhoE_x[:, 1:2]

    mass_flow_grad = tf.gradients(rho*u, x)[0]
    momentum_grad  = tf.gradients((rho*u*u + p), x)[0]
    energy_grad    = tf.gradients((rho*E + p)*u, x)[0]

    return [drho_t  + mass_flow_grad, 
            drhou_t + momentum_grad,
            drhoE_t + energy_grad,
    ]

geom = dde.geometry.Interval(0, 1)
timedomain = dde.geometry.TimeDomain(0, 2)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

# Initial conditions
def boundary(_, on_initial):
    return on_initial

def ic_rho(x):
    if x<0.5:
        rho = 1.4 * np.ones(x.shape)
    else:
        rho = 1.0 * np.ones(x.shape)

def ic_u(x):
    return 0.1 * np.ones(x.shape)

def ic_p(x):
    return 1.0 * np.ones(x.shape)

#ic1 = dde.IC(geomtime, lambda X: 1.4 * np.ones(X.shape) if (X.shape<0.5)  else 1.0 * np.ones(X.shape), boundary, component=0)
#ic2 = dde.IC(geomtime, lambda X: 0.1 * np.ones(X.shape),                                               boundary, component=1)
#ic3 = dde.IC(geomtime, lambda X: 1.0 * np.ones(X.shape),                                               boundary, component=2)

ic1 = dde.IC(geomtime, ic_rho, boundary, component=0)
ic2 = dde.IC(geomtime, ic_u,   boundary, component=1)
ic3 = dde.IC(geomtime, ic_p,   boundary, component=2)

# Boundary conditions
def boundary_l(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)

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

def bc_l_rho(x):
    return 1.4 * np.ones(x.shape)
def bc_l_u(x):
    return 0.1 * np.ones(x.shape)
def bc_l_p(x):
    return np.ones(x.shape)

def bc_r_rho(x):
    return np.ones(x.shape)
def bc_r_u(x):
    return 0.1 * np.ones(x.shape)
def bc_r_p(x):
    return np.ones(x.shape)

bc_l1 = dde.DirichletBC(geomtime, lambda X: 1.4 * np.ones(X.shape), boundary_l, component=0)
bc_l2 = dde.DirichletBC(geomtime, lambda X: 0.1 * np.ones(X.shape), boundary_l, component=1)
bc_l3 = dde.DirichletBC(geomtime, lambda X:         np.ones(X.shape), boundary_l, component=2)
bc_r1 = dde.DirichletBC(geomtime, lambda X:         np.ones(X.shape), boundary_r, component=0)
bc_r2 = dde.DirichletBC(geomtime, lambda X: 0.1  * np.ones(X.shape), boundary_r, component=1)
bc_r3 = dde.DirichletBC(geomtime, lambda X:         np.ones(X.shape), boundary_r, component=2)

#bc_l1 = dde.DirichletBC(geomtime, bc_l_rho, boundary_l, component=0)
#bc_l2 = dde.DirichletBC(geomtime, bc_l_u,   boundary_l, component=1)
#bc_l3 = dde.DirichletBC(geomtime, bc_l_p,   boundary_l, component=2)
#bc_r1 = dde.DirichletBC(geomtime, bc_r_rho, boundary_r, component=0)
#bc_r2 = dde.DirichletBC(geomtime, bc_r_u,   boundary_r, component=1)
#bc_r3 = dde.DirichletBC(geomtime, bc_r_p,   boundary_r, component=2)

data = dde.data.TimePDE(
    geomtime,
    Euler_system,
    [bc_l1, bc_l2, bc_l3, bc_r1, bc_r2, bc_r3, ic1, ic2, ic3],
    num_domain=400,
    num_boundary=40,
    num_initial=100,
    num_test=10000,
)

layer_size = [1] + [20] * 7 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)
model.compile("adam", lr=0.001, metrics=["l2 relative error"])

losshistory, train_state = model.train(epochs=10000)

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

checkpointer = dde.callbacks.ModelCheckpoint(
    "./model/model.ckpt", verbose=1, save_better_only=True
)
#movie = dde.callbacks.MovieDumper(
#    "model/movie", [-1], [1], period=100, save_spectrum=True, y_reference=solution
#)

# Plot PDE residue
x = geom.uniform_points(1000, True)
y = model.predict(x, operator=pde)
plt.figure()
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("PDE residue")
plt.show()

if name == "main":
main()

[lululxvi]

For example, we have a Dirichlet BC for the first output: u_0(x) = g(x). Then it is dde.DirichletBC(geomtime, g, ..., component=0). Assume the input X is N x d, then g should return an array of N x 1, but np.ones(X.shape) is N x d.

There are similar error in your IC.

Also, the following Python grammar is wrong,

def ic_rho(x):
    if x<0.5:
        rho = 1.4 * np.ones(x.shape)
    else:
        rho = 1.0 * np.ones(x.shape)

Here, input x is an numpy array. You cannot do x < 0.5.