Hi friends,
I’m trying to solve transient groundwater flow equation using PINNs. I have a rectangular domain 104 and the boundary conditions and initial conditions are presented below:
x-direction BCs: u(0, y, t) = 5 , u(10, y, t) = 1
y-direction BCs: du/dy = 0 @(x, 0, t) , du/dy = 0 @(x, 4, t)
Initial condition: u(x, y, 0) = -0.4x+5
I wrote the following code for this situation but I saw the water heads in the left and right boundaries are varying with time. Is it physically possible? I think it shouldn’t happen!
I’ve checked the code but I can’t find anything wrong.
Could you please help me find the problem with the code?
import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np
import matplotlib.pyplot as plt
Neural network model (simple feedforward network)
class Net(nn.Module):
def init(self):
super(Net, self).init()
self.fc1 = nn.Linear(3, 50)
self.fc2 = nn.Linear(50, 50)
self.fc3 = nn.Linear(50, 1)
def forward(self, x):
x = torch.tanh(self.fc1(x))
x = torch.tanh(self.fc2(x))
x = self.fc3(x)
return x
Physics-Informed Neural Network
alpha = 0.05 # hydraulic conductivity per Ss (0.00005/0.001)
mse = nn.MSELoss()
def flow_equation_loss(model, x, y, t):
# Compute prediction
u = model(torch.cat((x, y, t), dim=1))
# Gradients w.r.t. inputs (required for PDE)
u_t = torch.autograd.grad(u, t, torch.ones_like(u), create_graph=True, retain_graph=True)[0]
u_x = torch.autograd.grad(u, x, torch.ones_like(u), create_graph=True, retain_graph=True)[0]
u_xx = torch.autograd.grad(u_x, x, torch.ones_like(u_x), create_graph=True, retain_graph=True)[0]
u_y = torch.autograd.grad(u, y, torch.ones_like(u), create_graph=True, retain_graph=True)[0]
u_yy = torch.autograd.grad(u_y, y, torch.ones_like(u_y), create_graph=True, retain_graph=True)[0]
# Flow equation: u_t = alpha * (u_xx + u_yy)
pde_loss = mse(u_t, alpha * (u_xx + u_yy))
return pde_loss
def boundary_loss(model, x, y, t, u_b):
u_pred = model(torch.cat((x, y, t), dim=1))
return mse(u_pred, u_b)
def neumann_loss(model, x, y, t, target_flux, normal):
u = model(torch.cat((x, y, t), dim=1))
u_y = torch.autograd.grad(u, y, torch.ones_like(u), create_graph=True, retain_graph=True)[0]
flux_loss = mse(u_y * normal, target_flux)
return flux_loss
def train(model, x_interior, y_interior, t_interior,
x_boundary_left, y_boundary_left, t_boundary_left, u_boundary_left,
x_boundary_right, y_boundary_right, t_boundary_right, u_boundary_right,
x_neumann, y_neumann_top, y_neumann_bottom, t_neumann,
x_initial, y_initial, t_initial, u_initial,
optimizer, epochs=5000):
for epoch in range(epochs):
optimizer.zero_grad()
# PDE loss (interior points)
pde_loss = flow_equation_loss(model, x_interior, y_interior, t_interior)
# Boundary condition loss
boundary_loss_left = boundary_loss(model, x_boundary_left, y_boundary_left, t_boundary_left, u_boundary_left)
boundary_loss_right = boundary_loss(model, x_boundary_right, y_boundary_right, t_boundary_right, u_boundary_right)
# Neumann boundary condition loss for top and bottom (impermeable boundary, zero flux)
neumann_loss_top = neumann_loss(model, x_neumann, y_neumann_top, t_neumann, torch.zeros_like(x_neumann), normal=-1)
neumann_loss_bottom = neumann_loss(model, x_neumann, y_neumann_bottom, t_neumann, torch.zeros_like(x_neumann), normal=1)
# Initial condition loss
initial_loss = boundary_loss(model, x_initial, y_initial, t_initial, u_initial)
# Total loss
loss = pde_loss + 100*(boundary_loss_left + boundary_loss_right) + neumann_loss_top + neumann_loss_bottom + initial_loss
# Backpropagation
loss.backward(retain_graph=True)
optimizer.step()
if epoch % 500 == 0:
print(f"Epoch {epoch}, Loss: {loss.item()}")
print("Training completed!")
def predict(model, x, y, t):
return model(torch.cat((x, y, t), dim=1)).detach().numpy()
Training parameters
alpha = 0.05 # hydraulic conductivity per Ss (0.00005/0.001)
epochs = 5000
lr = 1e-3 # learning rate
Generate collocation points
N_interior = 5000 # number of points inside the domain
N_boundary = 1000 # number of points on the boundary
N_initial = 1000 # number of points for initial condition
Interior points (requires_grad=True for differentiation)
x_interior = 10torch.rand(N_interior, 1, requires_grad=True)
y_interior = 4torch.rand(N_interior, 1, requires_grad=True)
t_interior = torch.rand(N_interior, 1, requires_grad=True)
Boundary points for Dirichlet boundaries (left and right) with gradient tracking
x_boundary_left = 10torch.zeros(N_boundary // 2, 1, requires_grad=True)
y_boundary_left = 4torch.rand(N_boundary // 2, 1, requires_grad=True)
t_boundary_left = torch.rand(N_boundary // 2, 1, requires_grad=True)
x_boundary_right = 10torch.ones(N_boundary // 2, 1, requires_grad=True)
y_boundary_right = 4torch.rand(N_boundary // 2, 1, requires_grad=True)
t_boundary_right = torch.rand(N_boundary // 2, 1, requires_grad=True)
Boundary condition: u = 5 on the left boundary, u = 1 on the right boundary
u_boundary_left = torch.ones(N_boundary // 2, 1) * 5 # Left boundary at 5°C
u_boundary_right = torch.ones(N_boundary // 2, 1) * 1 # Right boundary at 1°C
Boundary points for Neumann boundaries (top and bottom) with gradient tracking
x_neumann = 10torch.rand(N_boundary // 2, 1, requires_grad=True)
y_neumann_top = 4torch.ones(N_boundary // 2, 1, requires_grad=True)
y_neumann_bottom = 4*torch.zeros(N_boundary // 2, 1, requires_grad=True)
t_neumann = torch.rand(N_boundary // 2, 1, requires_grad=True)
Initial condition: u(x, y, t=0) = -0.4*x + 5
x_initial = 10torch.rand(N_initial, 1, requires_grad=True)
y_initial = 4torch.rand(N_initial, 1, requires_grad=True)
t_initial = torch.zeros(N_initial, 1, requires_grad=True)
u_initial = -0.4 * x_initial + 5
Initialize the neural network and optimizer
model = Net()
optimizer = optim.Adam(model.parameters(), lr=lr)
Train the model
train(model, x_interior, y_interior, t_interior,
x_boundary_left, y_boundary_left, t_boundary_left, u_boundary_left,
x_boundary_right, y_boundary_right, t_boundary_right, u_boundary_right,
x_neumann, y_neumann_top, y_neumann_bottom, t_neumann,
x_initial, y_initial, t_initial, u_initial, optimizer, epochs)
Predict temperature distributions at 10 time steps
time_steps = np.linspace(0, 20, 10)
x_test = np.linspace(0, 10, 200) # finer mesh
y_test = np.linspace(0, 4, 200) # finer mesh
X, Y = np.meshgrid(x_test, y_test)
fig, axes = plt.subplots(2, 5, figsize=(20, 8)) # Plotting 10 figures
fig.subplots_adjust(hspace=0.3, wspace=0.3)
for i, t in enumerate(time_steps):
t_test = np.ones_like(X) * t
x_test_tensor = torch.tensor(X.reshape(-1, 1), dtype=torch.float32)
y_test_tensor = torch.tensor(Y.reshape(-1, 1), dtype=torch.float32)
t_test_tensor = torch.tensor(t_test.reshape(-1, 1), dtype=torch.float32)
u_pred = predict(model, x_test_tensor, y_test_tensor, t_test_tensor)
U_pred = u_pred.reshape(200, 200) # Use finer mesh for plotting
ax = axes.flat[i]
c = ax.contourf(X, Y, U_pred, 20, cmap='hot')
fig.colorbar(c, ax=ax)
ax.set_title(f"Time = {t:.2f}")
ax.set_xlabel("X")
ax.set_ylabel("Y")
plt.show()