multiprocessing working on progress

Battery.py

@@ -1,6 +1,11 @@
 from functools import reduce
 
+import logging
+import pickle
+import multiprocessing as mp
+from multiprocessing import Pool
 import numpy as np
+import numba as nb
 import pde
 import cantera as ct
 from scipy import optimize
@@ -17,6 +22,7 @@ class CombustionChamber:
         self.X0 = X0                            # Composition in mole fractions, Fuel + Air
         self.h0 = self.gas.enthalpy_mass        # inlet enthalpy
         self.hA = hA                            # HTC x Area
+        self.T1 = T0
         self.Twall0 = 1100 + 273.15
         self.Twall1 = 1100 + 273.15
         self.Area = 6.7 * 16.7
@@ -64,6 +70,19 @@ class CokeCharge:
     def end_baking (self, t):
         self.t_push = t
 
+brick_thickness = 0.14                            # m, 
+n_grid_brick = 16                                 # Number of Grid points inside
+wall_grid = pde.CartesianGrid([[0, brick_thickness]], n_grid_brick, periodic=False)
+wall_area = 6.7 * 16.7 # m^2 , Oven cross section area
+
+# op_grad2 = wall_grid.make_operator_no_bc('gradient_squared', backend='scipy')
+# op_grad = wall_grid.make_operator_no_bc('gradient', backend='scipy')
+# op_lap = wall_grid.make_operator_no_bc('laplace', backend='scipy')
+
+# op_info_grad2 = wall_grid._get_operator_info('gradient_squared')
+# op_info_grad = wall_grid._get_operator_info('gradient')
+# op_info_lap = wall_grid._get_operator_info('laplace')
+
 class CokeOvenBrickHeatEqn(pde.PDEBase):
     """Implementation of the Heat equation"""
 
@@ -90,20 +109,46 @@ class CokeOvenBrickHeatEqn(pde.PDEBase):
 
     def evolution_rate(self, state, t=0):
         """implement the python version of the evolution equation"""
-        state_lap = state.laplace(bc=self.bc)
-        state_grad = state.gradient(bc=self.bc)
-        
+        state_lap = state.laplace(bc=self.bc) # , backend="auto")
+        # state_grad = state.gradient(bc=self.bc, backend="scipy")
+        state_grad2 = state.gradient_squared(bc=self.bc) # , backend="auto")
+
+        '''
+        # out_cls_grad2 = state.get_class_by_rank(op_info_grad2.rank_out)
+        out_cls_grad = state.get_class_by_rank(op_info_grad.rank_out)
+        out_cls_lap = state.get_class_by_rank(op_info_lap.rank_out)
+        # state_grad2 = out_cls_grad2(state.grid, data="empty", dtype=state.dtype)
+        state_grad = out_cls_grad(state.grid, data="empty", dtype=state.dtype)
+        state_lap = out_cls_lap(state.grid, data="empty", dtype=state.dtype)
+        state.set_ghost_cells(self.bc)
+        # op_grad2(state._data_full, state_grad2.data)
+        op_grad(state._data_full, state_grad.data)
+        op_lap(state._data_full, state_lap.data)
+        '''
         k =  self.kCoef1  * state + self.kCoef0
         cp = self.cpCoef1 * state + self.cpCoef0
         
-        k_grad = self.kCoef1 * state_grad
+        state_grad_k_grad = self.kCoef1 * state_grad2 # state_grad.dot(state_grad)
         
-        return (state_grad.dot(k_grad) + k * state_lap) / self.rho / cp
+        return (state_grad_k_grad + k * state_lap) / cp / self.rho
 
-brick_thickness = 0.14                            # m, 
-n_grid_brick = 32                                 # Number of Grid points inside
-wall_grid = pde.CartesianGrid([[0, brick_thickness]], n_grid_brick, periodic=False)
-wall_area = 6.7 * 16.7 # m^2 , Oven cross section area
+    '''
+    def _make_pde_rhs_numba(self, state):
+        """implement the python version of the evolution equation"""
+        lap = state.grid.make_operator("laplace", bc=self.bc)
+        # grad = state.grid.make_operator("gradient", bc=self.bc)
+        grad2 = state.grid.make_operator("gradient_squared", bc=self.bc)
+        rho = self.rho
+        kCoef0 = self.kCoef0
+        kCoef1 = self.kCoef1
+        cpCoef0 = self.cpCoef0
+        cpCoef1 = self.cpCoef1
+        @nb.jit
+        def pde_rhs(data, t):
+            return (((kCoef1*grad2(data)) + (kCoef1*data + kCoef0)*lap(data)) / rho / (cpCoef1 * data + cpCoef0))
+        
+        return pde_rhs
+'''
 
 class RefractoryWall:
     def __init__ (self, T0):
@@ -125,7 +170,8 @@ class RefractoryWall:
 
     def solve (self, dt):
         # solution = self.eqn.solve (eqn, bc)
-        self.T_internal = self.eqn.solve(self.T_internal, t_range=dt, dt=1., tracker='consistency')
+        self.T_internal = self.eqn.solve(
+            self.T_internal, t_range=dt, dt=30., tracker='consistency', backend="numpy")
         self.T_chamber = self.T_internal.get_boundary_values(axis=0, upper=False, bc=self.eqn.bc)
 
     def heat_to_oven (self):
@@ -166,6 +212,10 @@ class OvenChamber:
         """ Update content with fresh coal is charged """
         self.content = coal_charge
 
+def wall_solve_wrapper(t_range, wall):
+    wall.solve(t_range)
+    return wall.T_internal, wall.T_chamber
+
 class Battery:
     
     def __init__ (self, name, size, heat_program, charge_program, burned_gas_state, hv):
@@ -185,6 +235,9 @@ class Battery:
         self.X0 = X0
         self.sequence_idx = 0                   # Integer, 0 ~ (size-1), progress index for oven sequence array
 
+        self.wall_t_history = []
+        self.gas_t_history = []
+
         self.hv = hv # Base unit heat J/kg
         self.normal_heat = self.heat_program.f(-1) # GJ / rev
 
@@ -220,7 +273,7 @@ class Battery:
 
         # 정상 상태 만들기: 모든 문에 n_cycle 회 장입
         n_cycle = 3             # 모든 문 장입 반복 횟수
-        period_over_dt = 11.    # period/dt, 장입 간격 / 초기화 time step 크기
+        period_over_dt = 6.    # period/dt, 장입 간격 / 초기화 time step 크기
         normal_period = self.charge_program.period(-1)  # 감산 전 장입 간격 (주기)
         dt =  normal_period / period_over_dt            # Simulation Time Step
         
@@ -229,6 +282,7 @@ class Battery:
         
         # initialization time loop
         for i in range(int(np.ceil(self.size * period_over_dt * n_cycle))):
+        # for i in range(3):
             """ Fill battety with normal charge rate """
             self.update(dt)     # Time adavancement
 
@@ -251,22 +305,22 @@ class Battery:
         # update chamber wall temperatures and mass flow rates
         # solve for equilibrium heat to walls
         for i_chamber, chmbr in enumerate(self.chambers):
-            try:
+            if i_chamber > 0:
                 wall_lower = self.walls_1[i_chamber-1]
-            except IndexError:
+            else:
                 wall_lower = None
 
-            try:
+            if i_chamber < self.size:
                 wall_upper = self.walls_0[i_chamber]
-            except IndexError:
+            else:
                 wall_upper = None
 
             chmbr.update_mdot(self.mdot(self.t))
             chmbr.update_Twall(
-                wall_lower.T_chamber if wall_lower else wall_upper.T_chamber, 
-                wall_upper.T_chamber if wall_upper else wall_lower.T_chamber, 
+                Twall0=(wall_lower.T_chamber if wall_lower else wall_upper.T_chamber), 
+                Twall1=(wall_upper.T_chamber if wall_upper else wall_lower.T_chamber), 
             )
-            print(f"t={self.t} - C{i_chamber} with {chmbr.Twall0} K and {chmbr.Twall1} K ")
+            print(f"t={self.t:6.2} : {chmbr.Twall0} K | Chamber {i_chamber} | {chmbr.Twall1} K ")
             chmbr.solve()
             Q1, Q2 = chmbr.heat() # W (J/s)
             if wall_lower: wall_lower.update_bc(Q=Q1)
@@ -282,13 +336,22 @@ class Battery:
             wall_lower.update_bc(T_oven=T_oven)
             wall_upper.update_bc(T_oven=T_oven)
 
-            wall_lower.solve(dt * 60 * 60)
-            wall_upper.solve(dt * 60 * 60)
+        with Pool(16) as pool:
+            wall_sln = pool.starmap(wall_solve_wrapper, [((dt*60*60), w) for w in self.walls_0+self.walls_1])
 
+            # wall_lower.solve(dt * 60 * 60) # convert hours to seconds
+            # wall_upper.solve(dt * 60 * 60) # convert hours to seconds
+        for ws, wall in zip(wall_sln, self.walls_0+self.walls_1):
+            T_internal, T_chamber = ws
+            wall.T_internal = T_internal
+            wall.T_chamber = T_chamber
+
+        '''
             ql = wall_lower.heat_to_oven()
             qu = wall_upper.heat_to_oven()
 
-            oven.bake(ql+qu)
+            oven.bake(ql+qu) 
+            '''
 
         # advance time oven brick
             # from chamber heat flux boundary condition
@@ -296,12 +359,6 @@ class Battery:
         
         # integrate heat to oven # 오븐 벽면 온도 우선 시간 함수로
 
-        '''dQ = self.dQ(dt) # array, dQ pairs of all oven taking from both walls
-        
-        for cc in self.processing:
-            cc.bake(dQ) # bake.(dQ[cc.idx_oven])
-            '''
-
     def push_and_charge (self, coke_charge):
         if len(self.processing) >= self.size:
             self.push(coke_charge.t_charge)
@@ -339,10 +396,6 @@ class Battery:
         # t_last + period 가 t, t + dt 사이에 들어오는 것 검사
         # t + dt 가 다음 추출/장입 시각 이후일 때 => 이번 time step 에 추출/장입을 실행해야함
         if self.t + dt >= period + self.t_last :
-            print(f"P/C within [ {self.t:7.3} , {self.t + dt:7.3} ].", 
-                  f"{self.t + dt - latest_coke_charge:7.3} since last P/C. ", 
-                  f"period = {self.charge_program.period(self.t):7.3}",)
-            
             # 마지막 장입 시각 + 장입 시간 간격 이 이번 time step 에 포함됨
             # 일정한 간격으로 장입 진행 중, 마지막 장입 시간 += 장입 간격
             if self.t < self.t_last + period:
@@ -357,10 +410,18 @@ class Battery:
             # oven = self.ovens[i_oven]
             fresh_coal = CokeCharge(self.t + dt, i_oven)
             self.push_and_charge(fresh_coal)
+
+            print(f"On {i_oven} P/C within [ {self.t:7.3} , {self.t + dt:7.3} ].", 
+                  f"{self.t + dt - latest_coke_charge:7.3} since last P/C. ", 
+                  f"period = {self.charge_program.period(self.t):7.3}",)
+            
         
         # 시뮬레이션 시간 업데이트
         self.t += dt
 
+        self.gas_t_history.append((self.t, [chmbr.T1 for chmbr in self.chambers]))
+        self.wall_t_history.append((self.t, [(wl.T_chamber, wl.T_internal.data, wl.T_oven, wu.T_oven, wu.T_internal.data, wu.T_chamber) for wl, wu in zip(self.walls_0, self.walls_1)]))
+
 def coke_oven_exhaust_stoichiometry (phi=1.0, return_unburned=False):
 
     # Define the oxidizer composition, here air with 21 mol-% O2 and 79 mol-% N2
@@ -487,4 +548,10 @@ if __name__ == "__main__":
     charging_plan = ChargeSchedule( 81, 9, 9, 1e-12, 24+13, 3, 1e-12 )
     n_doors = 66
 
-    bat3A = Battery("3A", n_doors, heating_plan, charging_plan, gas_in_state, hv)
+    bat3A = Battery("3A", n_doors, heating_plan, charging_plan, gas_in_state, hv)
+
+    with open('gas.history', 'wb') as gas_history_file:
+        pickle.dump(bat3A.gas_t_history, gas_history_file)
+
+    with open('wall.history', 'wb') as wall_history_file:
+        pickle.dump(bat3A.wall_t_history, wall_history_file)