1D shallow water

[yang1805524385]

Can DeepXDE simulate wave propagation under a submerged breakwater terrain controlled by a 1D shallow water equation PDE? I am looking forward to your reply.

[yang1805524385]

This is my code ，the results are incorrect; they seem to be unaffected by the terrain.
`device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
g = 9.81
xmin, xmax, tmax = 0., 30., 40.
input_dim = 2
output_dim = 2

def Slope(x):
if isinstance(x, np.ndarray): # 如果输入是 numpy 数组，转换为 torch.Tensor
x = torch.tensor(x, dtype=torch.float32)
y = torch.where(x < 7, 0.,
torch.where(x < 13, 0.05 * (x - 7),
torch.where(x < 15, 0.3,
torch.where(x < 18, 0.3 - 0.1 * (x - 15),
torch.where(x < 20, 0., 0.04 * (x - 20))))))
return y

Slope_x 函数：支持 PyTorch 张量

def Slope_x(x):
if isinstance(x, np.ndarray): # 如果输入是 numpy 数组，转换为 torch.Tensor
x = torch.tensor(x, dtype=torch.float32)
y = torch.where(x < 7, 0.,
torch.where(x < 13, 0.05,
torch.where(x < 15, 0.,
torch.where(x < 18, -0.1,
torch.where(x < 20, 0., 0.04)))))
return y

still water height

def sth(x):
if isinstance(x, np.ndarray): # 如果输入是 numpy 数组，转换为 torch.Tensor
x = torch.tensor(x, dtype=torch.float32)
y = torch.where(x < 7, 0.4,
torch.where(x < 13, 0.4-0.05 * (x - 7),
torch.where(x < 15, 0.1,
torch.where(x < 18, 0.1 + 0.1 * (x - 15),
torch.where(x < 20, 0.4, 0.4-0.04 * (x - 20))))))
return y

Define bottom topography

def get_bottom(bottom_name):
assert (isinstance(bottom_name, str))

bottom5 = lambda x: Slope(x)
bottom5_x = lambda x: Slope_x(x)

if bottom_name == 'b5':
    return {'b': bottom5, 'b_x': bottom5_x}
else:
    raise ValueError('Invalid bottom name, {} bottom not implemented'.format(bottom_name))

def _swe(bottom_name, g, source_term=True):
def swe(X, U):
bottom = get_bottom(bottom_name)
bx = bottom['b_x']
x, t = X[:, 0:1], X[:, 1:2]
h = U[:, 0:1]
u = U[:, 1:2]

    U1 = h
    U2 = h * u

    F1 = h * u
    F2 = h * u * u + 0.5 * g * h * h

    F1_x = dde.grad.jacobian(F1, X, i=0, j=0)
    F2_x = dde.grad.jacobian(F2, X, i=0, j=0)
    U1_t = dde.grad.jacobian(U1, X, i=0, j=1)
    U2_t = dde.grad.jacobian(U2, X, i=0, j=1)

    if source_term:
        h = 2 + torch.cos(x) * torch.cos(t)
        v = torch.sin(x) * torch.sin(t) / h
        S = torch.sin(x) * torch.cos(t) * (1 + v ** 2 - g * h) + 2 * v * torch.cos(x) * torch.sin(t) + g * h * bx(x)
    else:
        S = torch.zeros_like(x)
    b_x = bx(x)
    equaz_1 = U1_t + F1_x
    equaz_2 = U2_t + F2_x
    #+ g * h * b_x - S
    return [equaz_1, equaz_2]

return swe

def on_initial(_, on_initial):
return on_initial

def boundary(_, on_boundary):
return on_boundary

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

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

def func_ic_h(X):
if isinstance(X, np.ndarray): # 检查输入是否为 numpy.ndarray
X = torch.tensor(X, dtype=torch.float32)
x = X[:, 0:1]
h_ic = 0.4 + 0.1 * torch.sin(2 * torch.pi / 2.02 * x)
return h_ic

def _transform_output(U0):
def transform_output(x, U):
t = x[:, 1:2]
h = U[:, 0:1]
u = U[:, 1:2]
u_new = u * t + U0(x)
h_new = h
return torch.cat([h_new,u_new], axis=1)
return transform_output

def func_ic_u(x):
return 0.0

def func_bc_h1(X):
T = 2.02 # 周期 (秒)
A = 0.01 # 振幅 (cm)
if isinstance(X, np.ndarray): # 检查输入是否为 numpy.ndarray
X = torch.tensor(X, dtype=torch.float32)
t = X[:, 1:2]
#x_l = torch.ones_like(t) * xmin
h_lbc = 0.4 + A * torch.sin(2 * torch.pi / T * t)
return h_lbc

def func_bc_h2(X):
if isinstance(X, np.ndarray): # 检查输入是否为 numpy.ndarray
X = torch.tensor(X, dtype=torch.float32)
t = X[:, 1:2]
H2 = 0.43 + 0.01 * torch.sin(2 * torch.pi / 2.02 * t)
return 0.42

def func_bc_u1(X):
if isinstance(X, np.ndarray): # 检查输入是否为 numpy.ndarray
X = torch.tensor(X, dtype=torch.float32)
t = X[:, 1:2]
v = 0.0495torch.cos(3.1105t+90)
return v

def func_bc_h2(X):

if isinstance(X, np.ndarray): # 检查输入是否为 numpy.ndarray

X = torch.tensor(X, dtype=torch.float32)

H = ref_value[30, :]

print(H)

return H

model part

geom = dde.geometry.Interval(xmin, xmax)
timedomain = dde.geometry.TimeDomain(0.0, tmax)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

IC_h = dde.IC(geomtime, func_ic_h, on_initial, component=0)

IC_u = dde.IC(geomtime, func_ic_u, on_initial, component=1)

BC_h1 = dde.DirichletBC(geomtime, func_bc_h1, boundary_l, component=0)
BC_h2 = dde.DirichletBC(geomtime, func_bc_h2, boundary_r, component=0)

BC_h3 = dde.icbc.PointSetBC(12.5, ref_value[:, 37])

BC_u1 = dde.DirichletBC(geomtime, func_bc_u1, boundary_l, component=1)

BC = [IC_h, IC_u, BC_h1, BC_h2, BC_u1]
swe = _swe('b5', g, source_term=False)
data = dde.data.TimePDE(geomtime, swe, BC, num_domain=100, num_boundary=2, num_initial=100, num_test=100)
net = dde.maps.FNN(layer_sizes=[input_dim] + [60] * 5 + [output_dim], activation="tanh",
kernel_initializer="Glorot normal")
model = dde.Model(data, net)
resampler = dde.callbacks.PDEPointResampler(period=100)
model.compile('adam', lr=0.0005, decay = ("inverse time", 1000, 0.3), loss_weights=0.5)
transform_output = _transform_output(func_ic_u)
net.apply_output_transform(lambda x, U: transform_output(x, U))
losshistory, train_state = model.train(epochs=20_000)`