Comments (5)
lululxvi
commented on June 20, 2024
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.
from deepxde.
lkampoli
commented on June 20, 2024
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()
from deepxde.
lululxvi
commented on June 20, 2024
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
funcfor 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.
from deepxde.
lkampoli
commented on June 20, 2024
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()
from deepxde.
lululxvi
commented on June 20, 2024
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.
from deepxde.
Related Issues (20)
- Heat conduction equation with contact boundary and output-dependent parameters. HOT 3
- L-BFGS iteration records HOT 2
- How do I train pideepon with my own data? HOT 2
- How to solve inverse problem with parameters in the boundary condition? HOT 6
- Custom optimizer in deepxde HOT 1
- use a list of output layer in the PFNN HOT 2
- ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 4 and the array at index 1 has size 1
- No penetration Boundary condition HOT 2
- ZCS support for PINNs HOT 4
- Wave equation in a rectangular duct (waveguide)
- .
- L-BFGS HOT 1
- data files for deeponet examples HOT 1
- After I installed DeepxDE, guide the error JSON problem, and replace the deepxde version of the error error。json.decoder.JSONDecodeError: Expecting ',' delimiter: line 1 column 15 (char 14) HOT 3
- May I use the function like PointSetBc in physics-informed deeponet to feed some data for the model? HOT 1
- Specifying Dirichlet boundary condition using PointSetBC
- AttributeError: `dense` is not available with Keras 3. HOT 5
- How to set boundary condition which includes two PINN outputs
- Order of boundary points in on_boundary - Does it match the order in boundary_points argument of PointCloud? HOT 8
- For the L-BFGS optimizer, parameter "maxiter" should not be used to define the number of loops on data HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from deepxde.