Gen BCS working

GeneralizeBCS/Data/PINN/Tuner/results.txt

+... Tuner Results and Timings ...
+... Tuner timings ...
+Time elapsed: 33.408435106277466
+... Tuner results ...
+n_calls: 11
+Best error: 0.33532845973968506
+Best parameters: [3, 57, 'tanh']
+Best num_dense_layers: 3
+Best num_dense_nodes: 57
+Best activation: tanh

GeneralizeBCS/PINN2.py

+import deepxde as dde
+import numpy as np
+from deepxde.backend import tf
+import matplotlib.pyplot as plt
+
+from HeatCoefficients import *
+
+# Set random seed for reproducibility
+np.random.seed(42)
+tf.random.set_seed(42)
+
+# Global variables
+T_ICE_lb = 100.
+T_ICE_ub = 200.
+T_REF = 200.
+T_MELTING = 273.15
+Q_max = (T_ICE_ub-T_ICE_lb)/(T_MELTING-T_ICE_lb)
+delta_T = T_MELTING - T_ICE_lb  # ΔT
+beta = delta_T / T_ICE_lb  # β
+
+x1 = 1
+x0 = 0
+rho_ref = rho_physical(T_REF) # ρ_ref
+c_p_ref = c_p_physical(T_REF) # c_p_ref
+lambda_ref = lambda_physical(T_REF) # λ_ref
+L = x1 - x0
+t_ref = (rho_ref * c_p_ref * L**2 ) / lambda_ref # t_ref
+
+# Define the dimensionless PDE
+def pde(vars, theta):
+    # lhs
+    lhs = (
+        rho_star(theta, delta_T=delta_T, T_REF=T_REF, rho_ref=rho_ref) 
+        * c_p_star(theta, delta_T=delta_T, T_REF=T_REF, c_p_ref=c_p_ref) 
+        * dde.grad.jacobian(theta, vars, i=0, j=2)
+    )
+
+    # rhs
+    rhs = dde.grad.jacobian(
+        (lambda_star(theta, delta_T=delta_T, T_REF=T_REF, lambda_ref=lambda_ref) * dde.grad.jacobian(theta, vars, i=0, j=0) 
+        ),
+        vars, 
+        i=0, j=0
+    )
+    
+    # Dimensionless PDE Residual
+    return lhs - rhs
+
+def solve_PINN_dimensionless(t_end, layer, activation, initializer, num_domain=2_540, num_boundary=160, num_initial=160, solution=None, num_test=None, optimizer='adam', lr=1e-3, lr_min=1e-6, iterations=10_000, display_every=1_000, LBFGSB=True, savemodel='model.ckpt'):
+    # Define the domain
+    geom = dde.geometry.Rectangle([0, 0], [1, Q_max])
+    timedomain = dde.geometry.TimeDomain(0, t_end)
+    geomtime = dde.geometry.GeometryXTime(geom, timedomain)
+
+    # Initial condition: T(x,Q,0) = T_ICE
+    def initial_condition_full(X):
+        Q = X[:, 1:2]
+        return Q
+
+    ic = dde.icbc.IC(
+        geomtime,
+        initial_condition_full,
+        lambda X, on_initial: on_initial,
+        component=0
+    )
+
+    # Boundary Conditions: T(0, Q, tau) = T_ICE, T(1, Q, tau) = T_MELTING
+    def boundary_left_full(X, on_boundary):
+        return on_boundary and np.isclose(X[0], 0)
+
+    def boundary_right_full(X, on_boundary):
+        return on_boundary and np.isclose(X[0], 1)
+
+    def boundary_left_func(X):
+        Q = X[:, 1:2]
+        return Q
+
+    bc_left = dde.icbc.DirichletBC(
+        geomtime,
+        boundary_left_func,
+        boundary_left_full,
+        component=0
+    )
+
+    def boundary_right_func(X):
+        return np.ones_like(X[:, 0:1])
+
+    bc_right = dde.icbc.DirichletBC(
+        geomtime,
+        boundary_right_func,
+        boundary_right_full,
+        component=0
+    )
+
+    bcs = [ic, bc_left, bc_right]
+
+    # Define the data
+    data = dde.data.TimePDE(
+        geomtime, 
+        pde, 
+        bcs, 
+        num_domain=num_domain, 
+        num_boundary=num_boundary,
+        num_initial=num_initial,
+        solution=solution, 
+        num_test=num_test,
+    )
+
+    # Define NN
+    net = dde.maps.FNN(
+        layer, 
+        activation, 
+        initializer
+    )
+
+    # Apply hard bcs
+    net.apply_output_transform( # hard force initial condition
+        lambda x, y: x[:, 2:3] * (x[:, 0:1]) * y + x[:, 1:2]
+    )
+
+    # compile model
+    model = dde.Model(
+        data, 
+        net
+    )
+
+    # Save model while training
+    checkpointer = dde.callbacks.ModelCheckpoint(
+        savemodel, 
+        verbose=1,
+        save_better_only=True,
+    )
+
+    # Training loop
+    lr_current = lr
+    while lr_current >= lr_min:
+        print(f'current learnging rate: {lr_current}')
+        model.compile(
+            optimizer=optimizer, 
+            lr=lr_current,
+        )
+
+        losshistory, train_state = model.train(
+            iterations=iterations,
+            display_every=display_every,
+            callbacks=[checkpointer],
+        )
+        lr_current /= 10
+
+    # Use L-BFGS-B optimizer for better convergence
+    if LBFGSB:
+        model.compile("L-BFGS-B")
+        losshistory, train_state = model.train()
+
+    return model, losshistory, train_state
+
+
+
+# Plot the Loss History
+def plot_loss(losshistory):
+    dde.utils.plot_loss_history(losshistory)
+    plt.show()
+
+# Plot the Training State
+def plot_state(model, train_state):
+    dde.utils.plot_pde_residual(model, train_state)
+    plt.show()
+
+# Plot the Prediction
+def plot_prediction(model, x, Q, t):
+    # Create grid for x and t
+    X, T_grid = np.meshgrid(x, t)
+    X_flat = X.flatten()
+    T_flat = T_grid.flatten()
+    Q_flat = np.full_like(X_flat, Q)
+    X_input = np.vstack((X_flat, Q_flat, T_flat)).T
+    Y_pred = model.predict(X_input)
+    Y_pred = Y_pred.reshape(X.shape)
+
+    # Plot the prediction as a surface
+    fig = plt.figure(figsize=(10, 7))
+    ax = fig.add_subplot(111, projection='3d')
+    ax.plot_surface(X, T_grid, Y_pred)#, cmap='viridis')
+    ax.set_xlabel('x')
+    ax.set_ylabel('t')
+    ax.set_zlabel('T(x,t)')
+    ax.set_title(f'Heat Equation Solution for Q={Q}')
+    plt.show()
+
+def plot_prediction_1D(model, x, Q, t):
+    # Create grid for x and t
+    X, T_grid = np.meshgrid(x, t)
+    X_flat = X.flatten()
+    T_flat = T_grid.flatten()
+    Q_flat = np.full_like(X_flat, Q)
+    X_input = np.vstack((X_flat, Q_flat, T_flat)).T
+    Y_pred = model.predict(X_input)
+    Y_pred = Y_pred.reshape(X.shape)
+
+    # Plot the prediction as a surface
+    fig = plt.figure()
+    for i, time_ in enumerate(t):
+        plt.plot(x, Y_pred[i,:], label=f't={time_}')
+    plt.xlabel('x')
+    plt.ylabel('T(x,t)')
+    plt.legend()
+    plt.grid()
+    plt.show()
+
+def heat_flux_right(model, Q, t):
+    # Spatial derivative for PINN
+    def dThetadxi_PINN(x, y):
+        return dde.grad.jacobian(y, x, i=0, j=0)
+
+    # Create grid for x and t
+    x = np.linspace(0, 1, 100)
+    X, T_grid = np.meshgrid(x, t)
+    X_flat = X.flatten()
+    T_flat = T_grid.flatten()
+    Q_flat = np.full_like(X_flat, Q, dtype=float)
+    X_input = np.vstack((X_flat, Q_flat, T_flat)).T
+
+    dtheta_dxi_PINN = np.array(model.predict(X_input, operator=dThetadxi_PINN).reshape(len(t), len(x)))
+    Theta_PINN = np.array(model.predict(X_input).reshape(len(t), len(x)))
+    q_PINN = - dtheta_dxi_PINN[:, -1] * lambda_star(Theta_PINN[:, -1], delta_T, T_REF, lambda_ref)
+
+    return q_PINN 
+
+def heat_flux_right_physical(model, Q, t):
+    # Calculate heat flux
+    q_PINN = q_physical(heat_flux_right(model, Q, t), lambda_ref, delta_T, L)
+
+    return q_PINN 
+
+
+
+if __name__ == '__main__':
+    t_end = .1
+    layer = [3] + [20] * 3 + [1]
+    activation = "tanh"
+    initializer = "Glorot uniform"
+    model, losshistory, train_state = solve_PINN_dimensionless(
+        t_end=t_end,
+        layer=layer, activation=activation, initializer=initializer,
+        num_domain=2_540, num_boundary=160, num_initial=160,
+        num_test=2400,
+        optimizer='adam', lr=1e-3, lr_min=1e-3,
+        iterations=2_000, display_every=500,
+        savemodel='models/model.ckpt'
+    )
+
+    # Plot the loss history
+    plot_loss(losshistory)
+
+    # Plot the prediction
+    plot_prediction(
+        model=model, 
+        x=np.linspace(0,1,100), 
+        Q=.0, 
+        t=np.linspace(0, t_end, 100)
+    )
+
+    # Plot 1D Results
+    plot_prediction_1D(
+        model=model, 
+        x=np.linspace(0,1,100), 
+        Q=.0, 
+        t=np.linspace(0, t_end, 10)
+    )
+
+    # Calculate the heat flux at the right boundary and plot it for several Q values
+    # t = np.linspace(0, t_end_dimless, 100)
+    Q_values = np.linspace(0, Q_max, 10)
+    plt.figure()
+    t = np.linspace(0, t_end, 100)
+    for Q in Q_values:
+        q_PINN = heat_flux_right(
+            model=model, 
+            Q=Q, 
+            t=t
+        )
+        plt.plot(t, q_PINN, label=f'Q = {Q}')
+    plt.xlabel('t')
+    plt.ylabel('q(1, t)')
+    plt.legend()
+    plt.show()
+
+    # Calculate the heat flux (physical) at the right boundary and plot it for several Q values
+    plt.figure()
+    t = np.linspace(0, t_end, 100)
+    for Q in Q_values:
+        q_PINN = heat_flux_right_physical(
+            model=model, 
+            Q=Q, 
+            t=t
+        )
+        plt.plot(t, q_PINN, label=f'Q = {Q}')
+    plt.xlabel('t')
+    plt.ylabel('q(1, t)')
+    plt.legend()
+    plt.show()

GeneralizeBCS/TunerPINN.py

@@ -5,7 +5,7 @@ import deepxde as dde
 import tensorflow as tf
 import time
 
-from PINN import *
+from PINN2 import solve_PINN_dimensionless, plot_prediction, heat_flux_right
 from FEniCSx import solve_heatequation_dimensionless
 
 import skopt
@@ -14,6 +14,14 @@ from skopt.plots import plot_convergence, plot_objective
 from skopt.space import Real, Integer, Categorical, Integer
 from skopt.utils import use_named_args
 
+T_ICE_lb = 100.
+T_ICE_ub = 200.
+T_REF = 200.
+T_MELTING = 273.15
+Q_max = (T_ICE_ub-T_ICE_lb)/(T_MELTING-T_ICE_lb)
+delta_T = T_MELTING - T_ICE_lb  # ΔT
+beta = delta_T / T_ICE_lb  # β
+
 ####################################################################################################
 # Setup
 ####################################################################################################
@@ -21,35 +29,15 @@ from skopt.utils import use_named_args
 n_calls = 25
 
 # Physical properties
-T_ICE_lb = 100.
-T_ICe_ub = 200.
-T_MELTING = 273.15
 t_end = 1.
-# T_REF = 100.
-t_end_dimless = .1
-
-# Temperature Range
-delta_T = T_MELTING - T_ICE_lb  # ΔT = 100 K
 
-# Dimensionless Parameter
-beta = delta_T / T_ICE_lb        # β ≈ 0.3665
-
-# Reference Values at T_ICE
-rho_ref = rho_physical(T_REF)         # Reference density
-c_p_ref = c_p_physical(T_REF)         # Reference specific heat
-lambda_ref = lambda_physical(T_REF)   # Reference thermal conductivity
-
-# Reference Time Scale
-t_ref = (rho_ref * c_p_ref * L**2) / lambda_ref  # t_ref = (rho_ref * c_p_ref * L^2) / lambda_ref
-FO = (rho_ref * c_p_ref * L**2) / lambda_ref
-
-num_domain = 2540 # !10_000
-num_boundary = 80 # !100
-num_initial = 160 # !100
-num_test = 2540 # !None
+num_domain = 2_540 # !10_000
+num_boundary = 160 # !100
+num_initial = 240 # !100
+num_test = 2_000 # !None
 
 optimizer = 'adam'
-iterations = 5_000
+iterations_optimizer = 4_000
 display_every = 1_000
 learning_rate = 1e-3
 learning_rate_min = 1e-5
@@ -59,16 +47,12 @@ dim_num_dense_nodes = Integer(low=10, high=150, name='num_dense_nodes')
 dim_activation = Categorical(categories=['tanh', 'sigmoid', 'silu', 'swish'], name='activation')
 # dim_initializer = Categorical(categories=['Glorot uniform', 'Glorot uniform', 'He normal', 'He uniform'], name='initializer')
 initializer = 'Glorot uniform'
-# dim_learning_rate = Real(low=1e-4, high=1e-2, name='learning_rate', prior='log-uniform')
-# dim_T_REF = Integer(low=100, high=250, name='T_REF')
 
 dimensions = [
     dim_num_dense_layers,
     dim_num_dense_nodes,
     dim_activation,
     # dim_initializer,
-    # dim_learning_rate,
-    # dim_T_REF,
 ]
 
 default_parameters = [
@@ -76,8 +60,6 @@ default_parameters = [
     30, 
     'tanh', 
     # 'Glorot uniform',
-    # 1e-3, 
-    # 200,
 ]
 
 
@@ -96,36 +78,22 @@ def fitness(num_dense_layers, num_dense_nodes, activation):
     print('num_dense_nodes:', num_dense_nodes)
     print('activation:', activation)
     # print('initializer:', initializer)
-    # print('learning rate: {0:.1e}'.format(learning_rate))
     print()
 
     # Define ANN
     layer = [3] + [num_dense_nodes] * num_dense_layers + [1]
 
-    # model, losshistory, train_state = solve_PINN_dimensionless(
-    #     T_ICE, T_MELTING, t_end, T_REF, 
-    #     layer, activation, initializer, 
-    #     num_domain, num_boundary, num_initial, 
-    #     optimizer=optimizer, lr=learning_rate, 
-    #     iterations=iterations, display_every=display_every
-    # )
-    geomtime, bcs = geomXbcs(
-        T_ICE_lb, T_ICe_ub, 
-        T_MELTING, t_end_dimless
-    )
-    data, model = dataXmodel(
-        geomtime, bcs,
-        layer_size=layer, activation=activation, initializer=initializer,
-        num_domain=num_domain, num_boundary=num_boundary, 
-        num_initial=num_initial, num_test=num_test,
-        # optimizer=optimizer, lr=learning_rate
-    )
-
-    model, losshistory, train_state = train(
-        model, 
-        epochs=iterations,
-        display_every=display_every,
-        optimizer=optimizer, lr=learning_rate, lr_min=learning_rate_min
+    model, losshistory, train_state = solve_PINN_dimensionless(
+        t_end=t_end,
+        layer=layer, activation=activation, initializer=initializer,
+        num_domain=num_domain, 
+        num_boundary=num_boundary, num_initial=num_initial,
+        num_test=num_test,
+        optimizer=optimizer, 
+        lr=learning_rate, lr_min=learning_rate_min,
+        iterations=iterations_optimizer, display_every=display_every,
+        savemodel='models/PINN/TunerRuns/model.ckpt',
+        LBFGSB=True
     )
 
     error = np.array(losshistory.loss_train).sum(axis=1).ravel().min()
@@ -149,10 +117,10 @@ start_time_tuner = time.time()
 search_result = gp_minimize(
     func=fitness,
     dimensions=dimensions,
-    base_estimator='GP',#DUMMY', # ! 'GP', 'RF', 'ET', 'GBRT', 'DUMMY'
+    base_estimator='GP', # ! 'GP', 'RF', 'ET', 'GBRT', 'DUMMY'
     acq_func='EI', # Expected Improvement.
     n_calls=n_calls,
-    # n_random_starts=5, # ! set back to default
+    # n_random_starts=5, # Todo set back to default
     x0=default_parameters,
     random_state=42,
 )
@@ -160,7 +128,7 @@ end_time_tuner = time.time()
 print(f'Time elapsed: {end_time_tuner - start_time_tuner}')
 
 # Save results + timings
-with open('Data/PINNTunerTimings.txt', 'w') as f:
+with open('Data/PINN/Tuner/results.txt', 'w') as f:
     f.write('... Tuner Results and Timings ...\n')
     f.write(f'... Tuner timings ...\n')
     f.write(f'Time elapsed: {end_time_tuner - start_time_tuner}\n')
@@ -172,10 +140,9 @@ with open('Data/PINNTunerTimings.txt', 'w') as f:
     f.write(f'Best num_dense_nodes: {search_result.x[1]}\n')
     f.write(f'Best activation: {search_result.x[2]}\n')
     # f.write(f'Best initializer: {search_result.x[4]}\n')
-    # f.write(f'Best learning rate: {search_result.x[3]}\n')
 
 np.savez(
-    'Data/PINNTunerResults.npz', 
+    'Data/PINN/Tuner/results.npz', 
     time_elapsed=end_time_tuner - start_time_tuner, 
     n_calls=n_calls, best_error=search_result.fun, 
     best_parameters=search_result.x, 
@@ -183,10 +150,7 @@ np.savez(
     best_num_dense_nodes=search_result.x[1], 
     best_activation=search_result.x[2], 
     #best_initializer=search_result.x[4], 
-    # best_learning_rate=search_result.x[3],
     learning_rate=learning_rate, learning_rate_min=learning_rate_min,
-    T_ICE_lb=T_ICE_lb, T_ICE_ub=T_ICe_ub,
-    T_MELTING=T_MELTING
 )
 
 
@@ -195,10 +159,12 @@ print(search_result.x)
 
 # Plot convergence
 _ = plot_convergence(search_result)
-plt.savefig('Data/PINNTunerConvergence.png')
+plt.savefig('Data/PINN/Tuner/Convergence.png')
+plt.savefig('Data/PINN/Tuner/Convergence.pdf')
 plt.show()
 _ = plot_objective(search_result, show_points=True, size=3.8)
-plt.savefig('Data/PINNTunerObjective.png')
+plt.savefig('Data/PINN/Tuner/Objective.png')
+plt.savefig('Data/PINN/Tuner/Objective.pdf')
 plt.show()
 
 # Train best model
@@ -206,45 +172,51 @@ best_num_dense_layers=search_result.x[0]
 best_num_dense_nodes=search_result.x[1]
 best_activation=search_result.x[2]
 #best_initializer=search_result.x[4] 
-# best_learning_rate=search_result.x[3]
 
 best_layer = [3] + [best_num_dense_nodes] * best_num_dense_layers + [1]
-geomtime, bcs = geomXbcs(
-    T_ICE_lb, T_ICe_ub, 
-    T_MELTING, t_end_dimless
-)
-data, model = dataXmodel(
-    geomtime, bcs,
-    layer_size=best_layer, activation=best_activation, initializer=initializer,
-    num_domain=num_domain, num_boundary=num_boundary, 
-    num_initial=num_initial, num_test=num_test,
-    # optimizer=optimizer, lr=best_learning_rate
-)
 
-model, losshistory, train_state = train(
-    model, 
-    epochs=iterations,
-    display_every=display_every,
-    optimizer=optimizer, 
-    lr=learning_rate,
-    lr_min=learning_rate_min
-)
+model, losshistory, train_state = solve_PINN_dimensionless(
+        t_end=t_end,
+        layer=best_layer, activation=best_activation, 
+        initializer=initializer,
+        num_domain=num_domain, 
+        num_boundary=num_boundary, num_initial=num_initial,
+        num_test=num_test,
+        optimizer=optimizer, 
+        lr=learning_rate, lr_min=learning_rate_min,
+        iterations=iterations_optimizer, display_every=display_every,
+        savemodel='models/PINN/TunerBest/PINN.ckpt',
+        LBFGSB=True
+    )
 
 # Plot the Loss History
 dde.utils.plot_loss_history(losshistory)
-# dde.saveplot(
-#     losshistory, train_state, 
-#     issave=False, isplot=True
-# )
 plt.savefig('Data/PINNTunerLossHistory.png')
+plt.savefig('Data/PINNTunerLossHistory.pdf')
 plt.show()
 
 # Plot the prediction
-plot_prediction(
-    model=model, 
-    x=np.linspace(0,1,100), 
-    Q=.5, 
-    t=np.linspace(0, t_end_dimless, 100)
-)
-plt.savefig('Data/PINNTunerSolution_pcolormesh.png')
-plt.show()
+Q_values = np.linspace(0, Q_max, 10)
+for Q_ in Q_values:
+    # Plot the prediction
+    plot_prediction(
+        model=model, 
+        x=np.linspace(0,1,100), 
+        Q=Q_, 
+        t=np.linspace(0, t_end, 100)
+    )
+
+# heat flux at right boundary over time
+plt.figure()
+t = np.linspace(0, t_end, 100)
+for Q in Q_values:
+    q_PINN = heat_flux_right(
+        model=model, 
+        Q=Q, 
+        t=t
+    )
+    plt.plot(t, q_PINN, label=f'Q = {Q}')
+plt.xlabel('t')
+plt.ylabel('q(1, t)')
+plt.legend()
+plt.show()
\ No newline at end of file