aaa

Author: HTrun

examples/cylinder/3d_unsteady_continuous/load_lbm_data.py

@@ -159,9 +159,9 @@
 import copy
 from pyevtk.hl import pointsToVTK
 
-dir_ic = '/workspace/hesensen/PaddleScience_dev_3d/examples/cylinder/3d_unsteady_continuous/data/ic_data/'
+dir_ic = '/workspace/hesensen/PaddleScience_cqp/examples/cylinder/3d_unsteady_continuous/data/ic_data/'
 # dir_sp = '/home/aistudio/work/data/supervised_data/'
-dir_sp_mode1 = '/workspace/hesensen/PaddleScience_dev_3d/examples/cylinder/3d_unsteady_continuous/data/new_sp_data/'
+dir_sp_mode1 = '/workspace/hesensen/PaddleScience_cqp/examples/cylinder/3d_unsteady_continuous/data/new_sp_data/'
 
 #file_pattern = 'point_70.000000_9.000000_27.000000.dat'
 

examples/cylinder/3d_unsteady_continuous/test_3d_new.py

@@ -196,11 +196,14 @@
 # limitations under the License.
 
 import copy
+from abc import abstractmethod
+from typing import Union
 
 import numpy as np
 import paddle
 import paddle.distributed as dist
 import paddlescience as psci
+from paddle.optimizer import lr
 from pyevtk.hl import pointsToVTK
 
 import sample_boundary_training_data as sample_data
@@ -214,7 +217,7 @@
 # time arraep
 ic_t = 200000
 t_start = 200050
-t_end = 200250
+t_end = 200100
 t_step = 50
 time_num = int((t_end - t_start) / t_step) + 1
 time_tmp = np.linspace(t_start - ic_t, t_end - ic_t, time_num, endpoint=True)
@@ -229,76 +232,76 @@
 cylinder_wgt = 5.0
 top_wgt = 2.0
 bottom_wgt = 2.0
-eq_wgt= 10.0
+eq_wgt= 2.0
 ic_wgt = 5.0
 sup_wgt = 10.0
 
 # initial value
 txyz_uvwpe_ic = load_ic_data(ic_t)
-init_t = txyz_uvwpe_ic[:, 0]
-init_x = txyz_uvwpe_ic[:, 1]
-init_y = txyz_uvwpe_ic[:, 2]
-init_z = txyz_uvwpe_ic[:, 3]
-init_u = txyz_uvwpe_ic[:, 4]
-init_v = txyz_uvwpe_ic[:, 5]
-init_w = txyz_uvwpe_ic[:, 6]
-init_p = txyz_uvwpe_ic[:, 7]
+init_t = txyz_uvwpe_ic[:, 0]; print(f"init_t={init_t.shape} {init_t.mean().item():.10f}")
+init_x = txyz_uvwpe_ic[:, 1]; print(f"init_x={init_x.shape} {init_x.mean().item():.10f}")
+init_y = txyz_uvwpe_ic[:, 2]; print(f"init_y={init_y.shape} {init_y.mean().item():.10f}")
+init_z = txyz_uvwpe_ic[:, 3]; print(f"init_z={init_z.shape} {init_z.mean().item():.10f}")
+init_u = txyz_uvwpe_ic[:, 4]; print(f"init_u={init_u.shape} {init_u.mean().item():.10f}")
+init_v = txyz_uvwpe_ic[:, 5]; print(f"init_v={init_v.shape} {init_v.mean().item():.10f}")
+init_w = txyz_uvwpe_ic[:, 6]; print(f"init_w={init_w.shape} {init_w.mean().item():.10f}")
+init_p = txyz_uvwpe_ic[:, 7]; print(f"init_p={init_p.shape} {init_p.mean().item():.10f}")
 
 # num of supervised points
 n_sup = 2000
 
 # supervised data
 txyz_uvwpe_s = load_supervised_data(t_start, t_end, t_step, ic_t, n_sup)
-sup_t = txyz_uvwpe_s[:, 0]
-sup_x = txyz_uvwpe_s[:, 1]
-sup_y = txyz_uvwpe_s[:, 2]
-sup_z = txyz_uvwpe_s[:, 3]
-sup_u = txyz_uvwpe_s[:, 4]
-sup_v = txyz_uvwpe_s[:, 5]
-sup_w = txyz_uvwpe_s[:, 6]
-sup_p = txyz_uvwpe_s[:, 7]
+sup_t = txyz_uvwpe_s[:, 0]; print(f"sup_t={sup_t.shape} {sup_t.mean().item():.10f}")
+sup_x = txyz_uvwpe_s[:, 1]; print(f"sup_x={sup_x.shape} {sup_x.mean().item():.10f}")
+sup_y = txyz_uvwpe_s[:, 2]; print(f"sup_y={sup_y.shape} {sup_y.mean().item():.10f}")
+sup_z = txyz_uvwpe_s[:, 3]; print(f"sup_z={sup_z.shape} {sup_z.mean().item():.10f}")
+sup_u = txyz_uvwpe_s[:, 4]; print(f"sup_u={sup_u.shape} {sup_u.mean().item():.10f}")
+sup_v = txyz_uvwpe_s[:, 5]; print(f"sup_v={sup_v.shape} {sup_v.mean().item():.10f}")
+sup_w = txyz_uvwpe_s[:, 6]; print(f"sup_w={sup_w.shape} {sup_w.mean().item():.10f}")
+sup_p = txyz_uvwpe_s[:, 7]; print(f"sup_p={sup_p.shape} {sup_p.mean().item():.10f}")
 
 # num points to sample per GPU
-# num_points = 30000
-num_points = 15000
+num_points = 30000
+# num_points = 15000
 # discretize node by geo
 inlet_txyz, outlet_txyz, top_txyz, bottom_txyz, cylinder_txyz, interior_txyz = sample_data.sample_data(t_step=time_num, nr_points=num_points)
 
 # interior nodes discre
-i_t = interior_txyz[:, 0]
-i_x = interior_txyz[:, 1]
-i_y = interior_txyz[:, 2]
-i_z = interior_txyz[:, 3]
+i_t = interior_txyz[:, 0]; print(f"i_t={i_t.shape} {i_t.mean().item():.10f}")
+i_x = interior_txyz[:, 1]; print(f"i_x={i_x.shape} {i_x.mean().item():.10f}")
+i_y = interior_txyz[:, 2]; print(f"i_y={i_y.shape} {i_y.mean().item():.10f}")
+i_z = interior_txyz[:, 3]; print(f"i_z={i_z.shape} {i_z.mean().item():.10f}")
 
 # bc inlet nodes discre
-b_inlet_t = inlet_txyz[:, 0]
-b_inlet_x = inlet_txyz[:, 1]
-b_inlet_y = inlet_txyz[:, 2]
-b_inlet_z = inlet_txyz[:, 3]
+b_inlet_t = inlet_txyz[:, 0]; print(f"b_inlet_t={b_inlet_t.shape} {b_inlet_t.mean().item():.10f}")
+b_inlet_x = inlet_txyz[:, 1]; print(f"b_inlet_x={b_inlet_x.shape} {b_inlet_x.mean().item():.10f}")
+b_inlet_y = inlet_txyz[:, 2]; print(f"b_inlet_y={b_inlet_y.shape} {b_inlet_y.mean().item():.10f}")
+b_inlet_z = inlet_txyz[:, 3]; print(f"b_inlet_z={b_inlet_z.shape} {b_inlet_z.mean().item():.10f}")
 
 # bc outlet nodes discre
-b_outlet_t = outlet_txyz[:, 0]
-b_outlet_x = outlet_txyz[:, 1]
-b_outlet_y = outlet_txyz[:, 2]
-b_outlet_z = outlet_txyz[:, 3]
+b_outlet_t = outlet_txyz[:, 0]; print(f"b_outlet_t={b_outlet_t.shape} {b_outlet_t.mean().item():.10f}")
+b_outlet_x = outlet_txyz[:, 1]; print(f"b_outlet_x={b_outlet_x.shape} {b_outlet_x.mean().item():.10f}")
+b_outlet_y = outlet_txyz[:, 2]; print(f"b_outlet_y={b_outlet_y.shape} {b_outlet_y.mean().item():.10f}")
+b_outlet_z = outlet_txyz[:, 3]; print(f"b_outlet_z={b_outlet_z.shape} {b_outlet_z.mean().item():.10f}")
 
 # bc cylinder nodes discre
-b_cylinder_t = cylinder_txyz[:, 0]
-b_cylinder_x = cylinder_txyz[:, 1]
-b_cylinder_y = cylinder_txyz[:, 2]
-b_cylinder_z = cylinder_txyz[:, 3]
+b_cylinder_t = cylinder_txyz[:, 0]; print(f"b_cylinder_t={b_cylinder_t.shape} {b_cylinder_t.mean().item():.10f}")
+b_cylinder_x = cylinder_txyz[:, 1]; print(f"b_cylinder_x={b_cylinder_x.shape} {b_cylinder_x.mean().item():.10f}")
+b_cylinder_y = cylinder_txyz[:, 2]; print(f"b_cylinder_y={b_cylinder_y.shape} {b_cylinder_y.mean().item():.10f}")
+b_cylinder_z = cylinder_txyz[:, 3]; print(f"b_cylinder_z={b_cylinder_z.shape} {b_cylinder_z.mean().item():.10f}")
 
 # bc-top nodes discre
-b_top_t = top_txyz[:, 0] # value = [1, 2, 3, 4, 5]
-b_top_x = top_txyz[:, 1]
-b_top_y = top_txyz[:, 2]
-b_top_z = top_txyz[:, 3]
- 
+b_top_t = top_txyz[:, 0]; print(f"b_top_t={b_top_t.shape} {b_top_t.mean().item():.10f}") # value = [1, 2, 3, 4, 5]
+b_top_x = top_txyz[:, 1]; print(f"b_top_x={b_top_x.shape} {b_top_x.mean().item():.10f}")
+b_top_y = top_txyz[:, 2]; print(f"b_top_y={b_top_y.shape} {b_top_y.mean().item():.10f}")
+b_top_z = top_txyz[:, 3]; print(f"b_top_z={b_top_z.shape} {b_top_z.mean().item():.10f}")
+
 # bc-bottom nodes discre
-b_bottom_t = bottom_txyz[:, 0] # value = [1, 2, 3, 4, 5]
-b_bottom_x = bottom_txyz[:, 1]
-b_bottom_y = bottom_txyz[:, 2]
-b_bottom_z = bottom_txyz[:, 3]
+b_bottom_t = bottom_txyz[:, 0]; print(f"b_bottom_t={b_bottom_t.shape} {b_bottom_t.mean().item():.10f}") # value = [1, 2, 3, 4, 5]
+b_bottom_x = bottom_txyz[:, 1]; print(f"b_bottom_x={b_bottom_x.shape} {b_bottom_x.mean().item():.10f}")
+b_bottom_y = bottom_txyz[:, 2]; print(f"b_bottom_y={b_bottom_y.shape} {b_bottom_y.mean().item():.10f}")
+b_bottom_z = bottom_txyz[:, 3]; print(f"b_bottom_z={b_bottom_z.shape} {b_bottom_z.mean().item():.10f}")
 
 # bc & interior nodes for nn
 inputeq = np.stack((i_t, i_x, i_y, i_z), axis=1)
@@ -316,17 +319,17 @@
 pde = psci.pde.NavierStokes(nu=0.0205, rho=1.0, dim=3, time_dependent=True)
 
 # set bounday condition
-bc_inlet_u = psci.bc.Dirichlet("u", rhs=1)
-bc_inlet_v = psci.bc.Dirichlet("v", rhs=0)
-bc_inlet_w = psci.bc.Dirichlet("w", rhs=0)
-bc_cylinder_u = psci.bc.Dirichlet("u", rhs=0)
-bc_cylinder_v = psci.bc.Dirichlet("v", rhs=0)
-bc_cylinder_w = psci.bc.Dirichlet("w", rhs=0)
-bc_outlet_p = psci.bc.Dirichlet("p", rhs=0)
-bc_top_u = psci.bc.Dirichlet("u", rhs=0)
-bc_top_v = psci.bc.Dirichlet("v", rhs=0)
-bc_bottom_u = psci.bc.Dirichlet("u", rhs=0)
-bc_bottom_v = psci.bc.Dirichlet("v", rhs=0)
+bc_inlet_u = psci.bc.Dirichlet("u", rhs=1.0)
+bc_inlet_v = psci.bc.Dirichlet("v", rhs=0.0)
+bc_inlet_w = psci.bc.Dirichlet("w", rhs=0.0)
+bc_cylinder_u = psci.bc.Dirichlet("u", rhs=0.0)
+bc_cylinder_v = psci.bc.Dirichlet("v", rhs=0.0)
+bc_cylinder_w = psci.bc.Dirichlet("w", rhs=0.0)
+bc_outlet_p = psci.bc.Dirichlet("p", rhs=0.0)
+bc_top_u = psci.bc.Dirichlet("u", rhs=0.0)
+bc_top_v = psci.bc.Dirichlet("v", rhs=0.0)
+bc_bottom_u = psci.bc.Dirichlet("u", rhs=0.0)
+bc_bottom_v = psci.bc.Dirichlet("v", rhs=0.0)
 
 # add bounday and boundary condition
 pde.set_bc("inlet", bc_inlet_u, bc_inlet_v, bc_inlet_w)
@@ -345,8 +348,6 @@
 # Network
 net = psci.network.FCNet(
     num_ins=4, num_outs=4, num_layers=6, hidden_size=50, activation="tanh")
-# net.initialize("checkpoint/static_model_params_10000.pdparams")
-# net.initialize("checkpoint/dynamic_net_params_100000.pdparams")
 
 outeq = net(inputeq)
 outbc1 = net(inputbc1)
@@ -391,25 +392,161 @@
 algo = psci.algorithm.PINNs(net=net, loss=loss)
 
 # Optimizer
-opt = psci.optimizer.Adam(learning_rate=0.01, parameters=net.parameters())
+class LRBase(object):
+    """Base class for custom learning rates
+
+    Args:
+        epochs (int): total epoch(s)
+        step_each_epoch (int): number of iterations within an epoch
+        learning_rate (float): learning rate
+        warmup_epoch (int): number of warmup epoch(s)
+        warmup_start_lr (float): start learning rate within warmup
+        last_epoch (int): last epoch
+        by_epoch (bool): learning rate decays by epoch when by_epoch is True, else by iter
+        verbose (bool): If True, prints a message to stdout for each update. Defaults to False
+    """
+
+    def __init__(self,
+                 epochs: int,
+                 step_each_epoch: int,
+                 learning_rate: float,
+                 warmup_epoch: int,
+                 warmup_start_lr: float,
+                 last_epoch: int,
+                 by_epoch: bool,
+                 verbose: bool=False) -> None:
+        """Initialize and record the necessary parameters
+        """
+        super(LRBase, self).__init__()
+        if warmup_epoch >= epochs:
+            msg = f"When using warm up, the value of \"Global.epochs\" must be greater than value of \"Optimizer.lr.warmup_epoch\". The value of \"Optimizer.lr.warmup_epoch\" has been set to {epochs}."
+            print(msg)
+            warmup_epoch = epochs
+        self.epochs = epochs
+        self.step_each_epoch = step_each_epoch
+        self.learning_rate = learning_rate
+        self.warmup_epoch = warmup_epoch
+        self.warmup_steps = self.warmup_epoch if by_epoch else round(
+            self.warmup_epoch * self.step_each_epoch)
+        self.warmup_start_lr = warmup_start_lr
+        self.last_epoch = last_epoch
+        self.by_epoch = by_epoch
+        self.verbose = verbose
+
+    @abstractmethod
+    def __call__(self, *kargs, **kwargs) -> lr.LRScheduler:
+        """generate an learning rate scheduler
+
+        Returns:
+            lr.LinearWarmup: learning rate scheduler
+        """
+        pass
+
+    def linear_warmup(
+            self,
+            learning_rate: Union[float, lr.LRScheduler]) -> lr.LinearWarmup:
+        """Add an Linear Warmup before learning_rate
+
+        Args:
+            learning_rate (Union[float, lr.LRScheduler]): original learning rate without warmup
+
+        Returns:
+            lr.LinearWarmup: learning rate scheduler with warmup
+        """
+        warmup_lr = lr.LinearWarmup(
+            learning_rate=learning_rate,
+            warmup_steps=self.warmup_steps,
+            start_lr=self.warmup_start_lr,
+            end_lr=self.learning_rate,
+            last_epoch=self.last_epoch,
+            verbose=self.verbose)
+        return warmup_lr
+
+
+class Constant(lr.LRScheduler):
+    """Constant learning rate Class implementation
+
+    Args:
+        learning_rate (float): The initial learning rate
+        last_epoch (int, optional): The index of last epoch. Default: -1.
+    """
+
+    def __init__(self, learning_rate, last_epoch=-1, **kwargs):
+        self.learning_rate = learning_rate
+        self.last_epoch = last_epoch
+        super(Constant, self).__init__()
+
+    def get_lr(self) -> float:
+        """always return the same learning rate
+        """
+        return self.learning_rate
+
+
+class Cosine(LRBase):
+    """Cosine learning rate decay
+
+    ``lr = 0.05 * (math.cos(epoch * (math.pi / epochs)) + 1)``
+
+    Args:
+        epochs (int): total epoch(s)
+        step_each_epoch (int): number of iterations within an epoch
+        learning_rate (float): learning rate
+        eta_min (float, optional): Minimum learning rate. Defaults to 0.0.
+        warmup_epoch (int, optional): The epoch numbers for LinearWarmup. Defaults to 0.
+        warmup_start_lr (float, optional): start learning rate within warmup. Defaults to 0.0.
+        last_epoch (int, optional): last epoch. Defaults to -1.
+        by_epoch (bool, optional): learning rate decays by epoch when by_epoch is True, else by iter. Defaults to False.
+    """
+
+    def __init__(self,
+                 epochs,
+                 step_each_epoch,
+                 learning_rate,
+                 eta_min=0.0,
+                 warmup_epoch=0,
+                 warmup_start_lr=0.0,
+                 last_epoch=-1,
+                 by_epoch=False,
+                 **kwargs):
+        super(Cosine, self).__init__(epochs, step_each_epoch, learning_rate,
+                                     warmup_epoch, warmup_start_lr, last_epoch,
+                                     by_epoch)
+        self.T_max = (self.epochs - self.warmup_epoch) * self.step_each_epoch
+        self.eta_min = eta_min
+        if self.by_epoch:
+            self.T_max = self.epochs - self.warmup_epoch
+
+    def __call__(self):
+        learning_rate = lr.CosineAnnealingDecay(
+            learning_rate=self.learning_rate,
+            T_max=self.T_max,
+            eta_min=self.eta_min,
+            last_epoch=self.last_epoch) if self.T_max > 0 else Constant(
+                self.learning_rate)
+
+        if self.warmup_steps > 0:
+            learning_rate = self.linear_warmup(learning_rate)
+
+        setattr(learning_rate, "by_epoch", self.by_epoch)
+        return learning_rate
+
+num_epoch = 100000
+_lr = Cosine(num_epoch, 1, 0.001, warmup_epoch=5000, by_epoch=True)()
+# _lr = 0.001
+opt = psci.optimizer.Adam(learning_rate=_lr, parameters=net.parameters())
 
 # Solver
 solver = psci.solver.Solver(pde=pde, algo=algo, opt=opt)
 
 # Solve
-solution = solver.solve(num_epoch=100000)
-# solution = solver.predict()
-# print(type(solution), len(solution), solution[0].shape)
-# exit()
-
-# solution = [
-#     np.stack([init_u, init_v, init_w, init_p], axis=1)
-# ]
+solution = solver.solve(num_epoch=num_epoch)
+
+# print shape of every subset points' output
 for idx, si in enumerate(solution):
     print(f"solution[{idx}].shape = {si.shape}")
 
+# only coord at start time is needed
 n = int(i_x.shape[0] / len(time_array))
-# n = 1
 i_x = i_x.astype("float32")
 i_y = i_y.astype("float32")
 i_z = i_z.astype("float32")
@@ -420,6 +557,5 @@
 i_z = i_z * 320
 
 cord = np.stack((i_x[0:n], i_y[0:n], i_z[0:n]), axis=1)
-# cord = np.stack((i_x, i_y, i_z), axis=1)
 # psci.visu.__save_vtk_raw(cordinate=cord, data=solution[0][-n::])
-psci.visu.save_vtk_cord(filename="./vtk/output_2023_1_12", time_array=time_array, cord=cord, data=solution)
+psci.visu.save_vtk_cord(filename="./vtk/output_2023_1_13_new", time_array=time_array, cord=cord, data=solution)