# -*- coding: utf-8 -*-
"""
Created on Tue Apr  2 15:14:40 2019

@author: ddesantis
"""

from fipy import *
import numpy as np
import matplotlib.pyplot as plt

#%% - Mesh

nx = 50
L = 0.061 #m,

L1 = 0.035              # Length of air channel
t_w = 1./1000.          # m, wall thickness
L2 = L1+t_w             # m, interface between PCM and wall
L3 = L                  # m, Total Length
mesh = Grid1D(Lx=L,nx=nx)
#mesh = CylindricalGrid1D(nr=nx, dr=dx)

X = mesh.faceCenters
x=mesh.cellCenters[0]
#%% - Sweeps
dt = 0.45                               #s,time step for sweeps
steps = 150
#%% Universal Constants
MW = 2.02        #g H2 mol^-1
R = 8.314e-3    # kj mol^-1 K^-1

#%% - air Properties
T_in = 600.      # K
P_in = 1.       # bar
v = 0.5         # m s^-2
k_air = 44.41/1000./1000.       # kJ s^-1 m^-1 K^-1 (air at ~600K)
Cp_air = 1.051  # kJ kg ^-1 K^-1
rho_air = 0.6   # kg m^-3, At 600K

Re = 1150.   # Reynolds Number, dimensionless
Pr = 0.727  # Prandlt Number, dimensionless
Nu = 1.86   # Nusselt Number, dimensionless

#%% - Wall properties
#I'm assuming the wall is steel
rho_w = 8050.     # kg m^-3
Cp_w = 0.5       # kJ K^-1 kg^-1
k_w = 16.3/1000.  # kJ s^-1 m^-1 K^-1

T_w_initial = (T_in+293)/2

#%% PCM Properties
T_init = 293.    #K
T_melt = 300.    #K
H_fus = 206.     #kJ kg^-1

k_s = 0.18      #W m^-1 K^-1
k_l = 0.19      #W m^-1 K^-1
#k_PCM = FaceVariable(mesh=mesh,value = k_s)
k_PCM = k_s

Cp_s = 1.8      #kJ kg^-1 K^-1
Cp_l = 2.4      #kJ kg^-1 K^-1
#Cp_PCM = FaceVariable(mesh=mesh,value = Cp_s)
Cp_PCM = Cp_s

# I assume they want solid and liquid presented the same way again?
rho_s = 789.     # kg m^-3
rho_l = 750.     # kg m^-3

#%% Domains

air = x <= L1
wall = (x > L1) & (x <= L2)
PCM = x > L2

#%% Equation Variables

T = CellVariable(name='Temp', mesh=mesh, hasOld=True)
T.setValue(T_in, where=air)
T.setValue(T_w_initial, where=wall)
T.setValue(T_init, where=PCM)

rho = CellVariable(name='rho', mesh=mesh)
rho.setValue(rho_air, where=air)
rho.setValue(rho_w, where=wall)
rho.setValue(rho_s, where=PCM)

Cp = CellVariable(name='Cp', mesh=mesh)
Cp.setValue(Cp_air, where=air)
Cp.setValue(Cp_w, where=wall)
Cp.setValue(Cp_s, where=PCM)

k = CellVariable(name='k', mesh=mesh)
k.setValue(k_air, where=air)
k.setValue(k_w, where=wall)
k.setValue(k_s, where=PCM)

#%% Boundary Conditions

T.constrain(T_in, mesh.facesLeft)

#%% - Equations

eq = TransientTerm(coeff=rho * Cp, var=T) == DiffusionTerm(coeff=k, var=T)

#%% - Sweeps
viewer = Viewer(vars=T,datamin = 50.,datamax = 700.)
viewer.plot()

for step in range(steps):
    T.updateOld()
    eq.solve(dt=dt)

    viewer.plot()
    percent_complete = float(step)/float(steps)
    print('percent complete: {:0.0f}%'.format(100.*percent_complete))

    print('step {} of {}'.format(step, steps))

    print('elapsed fill time: {:0.2f}s'.format(step*dt))