Pull Request: Implement PINN Solver for Timoshenko Beam Theory with Two-Stage Optimization

examples/pinn_forward/Euler_beam_noisy_measurements.py

@@ -0,0 +1,158 @@
+"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch, paddle"""
+import deepxde as dde
+import numpy as np
+
+
+def ddy(x, y):
+    return dde.grad.hessian(y, x)
+
+
+def dddy(x, y):
+    return dde.grad.jacobian(ddy(x, y), x)
+
+
+def pde(x, y):
+    dy_xx = ddy(x, y)
+    dy_xxxx = dde.grad.hessian(dy_xx, x)
+    return dy_xxxx + 1
+
+
+def boundary_l(x, on_boundary):
+    return on_boundary and dde.utils.isclose(x[0], 0)
+
+
+def boundary_r(x, on_boundary):
+    return on_boundary and dde.utils.isclose(x[0], 1)
+
+
+# Closed-form solution
+def func(x):
+    return -(x**4) / 24 + x**3 / 6 - x**2 / 4
+
+
+# ---------------------------------------------------------
+# Synthetic FEM / experimental measurements
+# ---------------------------------------------------------
+
+np.random.seed(42)
+
+num_measurements = 10
+
+# Random measurement locations
+X_measurement = np.random.uniform(
+    0, 1, (num_measurements, 1)
+)
+
+# Exact displacement at measurement locations
+Y_exact = func(X_measurement)
+
+# Simulated numerical/measurement error
+noise_level = 0.005  # 0.5%
+
+noise = np.random.normal(
+    loc=0.0,
+    scale=noise_level,
+    size=Y_exact.shape
+)
+
+# FEM-like noisy measurements
+Y_measurement = Y_exact * (1.0 + noise)
+
+
+# Treat measurements as observed displacement data
+measurement_bc = dde.icbc.PointSetBC(
+    X_measurement,
+    Y_measurement,
+    component=0,
+)
+
+
+# ---------------------------------------------------------
+# Geometry and boundary conditions
+# ---------------------------------------------------------
+
+geom = dde.geometry.Interval(0, 1)
+
+bc1 = dde.icbc.DirichletBC(
+    geom,
+    lambda x: 0,
+    boundary_l
+)
+
+bc2 = dde.icbc.NeumannBC(
+    geom,
+    lambda x: 0,
+    boundary_l
+)
+
+bc3 = dde.icbc.OperatorBC(
+    geom,
+    lambda x, y, _: ddy(x, y),
+    boundary_r
+)
+
+bc4 = dde.icbc.OperatorBC(
+    geom,
+    lambda x, y, _: dddy(x, y),
+    boundary_r
+)
+
+
+# ---------------------------------------------------------
+# PINN data
+# ---------------------------------------------------------
+
+data = dde.data.PDE(
+    geom,
+    pde,
+    [
+        bc1,
+        bc2,
+        bc3,
+        bc4,
+        measurement_bc,
+    ],
+    num_domain=10,
+    num_boundary=2,
+    solution=func,
+    num_test=100,
+)
+
+
+# ---------------------------------------------------------
+# Neural network
+# ---------------------------------------------------------
+
+layer_size = [1] + [20] * 3 + [1]
+activation = "tanh"
+initializer = "Glorot uniform"
+
+net = dde.nn.FNN(
+    layer_size,
+    activation,
+    initializer
+)
+
+
+# ---------------------------------------------------------
+# Training
+# ---------------------------------------------------------
+
+model = dde.Model(data, net)
+
+model.compile(
+    "adam",
+    lr=0.001,
+    metrics=["l2 relative error"]
+)
+
+losshistory, train_state = model.train(
+    iterations=10000
+)
+
+dde.saveplot(
+    losshistory,
+    train_state,
+    issave=True,
+    isplot=True
+)

examples/pinn_forward/Timoshenko_beam.py

@@ -0,0 +1,97 @@
+"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch, paddle"""
+import deepxde as dde
+import matplotlib.pyplot as plt
+import numpy as np
+
+EI = 1.0
+GA_s = 5.0
+q = -1.0
+
+def pde(x, y):
+    # y: [w, phi]
+    dw_x = dde.grad.jacobian(y, x, i=0, j=0)
+    dphi_x = dde.grad.jacobian(y, x, i=1, j=0)
+
+    dw_xx = dde.grad.hessian(y, x, component=0)
+    dphi_xx = dde.grad.hessian(y, x, component=1)
+
+    res_w = GA_s * (dw_xx - dphi_x) + q
+    res_phi = EI * dphi_xx + GA_s * (dw_x - y[:, 1:2])
+
+    return [res_w, res_phi]
+
+def boundary_l(x, on_boundary):
+    return on_boundary and dde.utils.isclose(x[0], 0)
+
+def boundary_r(x, on_boundary):
+    return on_boundary and dde.utils.isclose(x[0], 1)
+
+# Boundary conditions: Free end M=0 and V=0
+def bc_moment_free(x, y, _):
+    return dde.grad.jacobian(y, x, i=1, j=0)
+
+def bc_shear_free(x, y, _):
+    dw_x = dde.grad.jacobian(y, x, i=0, j=0)
+    return dw_x - y[:, 1:2]
+
+# Analytical solution
+def func(x):
+    w_b = -(x**4) / 24.0 + (x**3) / 6.0 - (x**2) / 4.0
+    w_s = (1.0 / GA_s) * ((x**2) / 2.0 - x)
+    w_total = w_b + w_s
+    phi = -(x**3) / 6.0 + (x**2) / 2.0 - x / 2.0
+    return np.hstack((w_total, phi))
+
+
+geom = dde.geometry.Interval(0, 1)
+
+bc_w0 = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_l, component=0)
+bc_phi0 = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_l, component=1)
+bc_ML = dde.icbc.OperatorBC(geom, bc_moment_free, boundary_r)
+bc_VL = dde.icbc.OperatorBC(geom, bc_shear_free, boundary_r)
+
+data = dde.data.PDE(
+    geom,
+    pde,
+    [bc_w0, bc_phi0, bc_ML, bc_VL],
+    num_domain=60,
+    num_boundary=2,
+    solution=func,
+    num_test=100,
+)
+
+net = dde.nn.FNN([1] + [30] * 3 + [2], "tanh", "Glorot uniform")
+model = dde.Model(data, net)
+
+# Weights ranking: [res_w, res_phi, bc_w0, bc_phi0, bc_ML, bc_VL]
+loss_weights = [1.0, 1.0, 10.0, 10.0, 5.0, 5.0]
+
+# 1. stage: Adam convergence
+#model.compile("adam", lr=0.001, loss_weights=loss_weights ,metrics=["l2 relative error"])
+model.compile("adam", lr=0.001, metrics=["l2 relative error"])
+model.train(iterations=6000)
+
+# 2. stage: L-BFGS fine tuning
+#model.compile("L-BFGS", loss_weights=loss_weights, metrics=["l2 relative error"])
+model.compile("L-BFGS", metrics=["l2 relative error"])
+losshistory, train_state = model.train()
+
+# Visualization
+x_test = np.linspace(0, 1, 200).reshape(-1, 1)
+y_pred = model.predict(x_test)
+y_true = func(x_test)
+
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
+ax1.plot(x_test, y_true[:, 0], "k--", label="Analytical (Exact)")
+ax1.plot(x_test, y_pred[:, 0], "r-", label="PINN prediction")
+ax1.set_title("Transverse displacement $w(x)$")
+ax1.grid(True, linestyle=":", alpha=0.6)
+ax1.legend()
+
+ax2.plot(x_test, y_true[:, 1], "k--", label="Analytical (Exact)")
+ax2.plot(x_test, y_pred[:, 1], "b-", label="PINN prediction")
+ax2.set_title(r"Section rotation $\varphi(x)$")
+ax2.grid(True, linestyle=":", alpha=0.6)
+ax2.legend()
+plt.tight_layout()
+plt.show()