diff --git a/examples/pinn_forward/Euler_beam_noisy_measurements.py b/examples/pinn_forward/Euler_beam_noisy_measurements.py new file mode 100644 index 000000000..38c943e12 --- /dev/null +++ b/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 +) \ No newline at end of file diff --git a/examples/pinn_forward/Timoshenko_beam.py b/examples/pinn_forward/Timoshenko_beam.py new file mode 100644 index 000000000..4f2044763 --- /dev/null +++ b/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()