# Reverse-Engineering Thermal Diffusivity in Mulch

The insulation considered for over-summer snow storage uses a thick layer of mulch. From the perspective of heat transfer, mulch is a complex material which can neither be described as homogeneous nor as porous. The local void volume fraction varies significantly and can be of the order of the mulch volume fraction. The mulch chips are anisotropic due to the natural alignment of wood fibers and porous. In the application of interest here, they hold significant humidity.

Not surprisingly, there is no known correlation for the apparent thermal conductivity of mulch. This notebook investigates the possibility of reverse-engineering this critical property from the temperature data collected at the Craftsbury site.

## Equations

We solve the unsteady heat equation over the insulation, assuming 1D heat transfer. The mulch is to be modeled by an apparent thermal diffusivity $\alpha$, which lumps together conduction in the wood chips, possible convection in the air cavities and evaporation. The resulting equation is

$$
\frac{\partial T}{\partial t}=\frac{\partial}{\partial x}\alpha\frac{\partial T}{\partial x}+q''
$$

where $\alpha$ may be a function of multiple variables including local temperature, depth, humidity, ...

$$
q'' = -h\Delta T
$$

$$
q'' = -h(T_f,\Delta T)\Delta T
$$

## Approach

We seek to solve this equation using the temperature at the top (underneath the reflective blanket) and the bottom (on the surface of the snow) of the insulation. The temperature signal at $10cm$ from the snow serves as the target our simulation should predict to a reasonable uncertainty. 

The above heat equation is discretized on a non-uniform 1D grid using second order finite differences and the Crank Nicolson time advancement method:

$$
\frac{T_i^{(n+1)}-T_i^{(n)}}{\Delta T}=\frac{1}{2}\left(\left.\frac{\delta}{\delta x}\alpha\frac{\delta T}{\delta x}\right\vert_i^{(n+1)}+\left.\frac{\delta}{\delta x}\alpha\frac{\delta T}{\delta x}\right\vert_i^{(n)}\right)+{q''}^{(n)}_i
$$

$$
\frac{T_i^{(n+1)}-T_i^{(n)}}{\Delta T}=\frac{1}{2}\left(a_iT_{i-1}^{(n+1)}+b_iT_{i}^{(n+1)}+c_iT_{i+1}^{(n+1)}+a_iT_{i-1}^{(n)}+b_iT_{i}^{(n)}+c_iT_{i+1}^{(n)}\right)+{q''}^{(n)}_i
$$

which is recast into the following linear system:

$$
-a_i\frac{\Delta t}{2}T_{i-1}^{(n+1)}+\left(1 -b_i\frac{\Delta t}{2}\right)T_{i}^{(n+1)}-c_i\frac{\Delta t}{2}T_{i+1}^{(n+1)}=a_i\frac{\Delta t}{2}T_{i-1}^{(n)}+\left(1+b_i\frac{\Delta t}{2}T_{i}^{(n)}\right)+c_i\frac{\Delta t}{2}T_{i+1}^{(n)}+\frac{\Delta t}{2}{q''}^{(n)}_i
$$

The system is solved as an inversion matrix $A.x=d$ where $A$ is a tridiagonal matrix of lower, main and upper diagonals $a_i$, $b_i$, $c_i$, respectively, $x$ is the vector of temperature nodes at time $n+1$ $x=[T_0^{(n+1)},...,T_i^{(n+1)},...,T_{nx-1}^{(n+1)}]^T$ and $d$ is a vector made of quantities calculated at time $n$.

To reverse engineer the diffusivity of mulch, we first assume that the diffusivities are constant and we can integrate the above equation in time using data collected over last summer.



### Python set-up and useful functions


```python
%matplotlib inline 
# plots graphs within the notebook

from IPython.display import display,Image, Latex
from sympy.interactive import printing
printing.init_printing(use_latex='mathjax')
from IPython.display import clear_output

import time

import numericaltools as numtools


import matplotlib.pyplot as plt
import numpy as np
import math
import scipy.constants as sc
import h5py


import sympy as sym




class PDF(object):
  def __init__(self, pdf, size=(200,200)):
    self.pdf = pdf
    self.size = size

  def _repr_html_(self):
    return ''.format(self.pdf, self.size)

  def _repr_latex_(self):
    return r'\includegraphics[width=1.0\textwidth]{{{0}}}'.format(self.pdf)

class ListTable(list):
    """ Overridden list class which takes a 2-dimensional list of 
        the form [[1,2,3],[4,5,6]], and renders an HTML Table in 
        IPython Notebook. """
    
    def _repr_html_(self):
        html = ["<table>"]
        for row in self:
            html.append("<tr>")
            
            for col in row:
                html.append("<td>{0}</td>".format(col))
            
            html.append("</tr>")
        html.append("</table>")
        return ''.join(html)
    
font = {'family' : 'serif',
        #'color'  : 'black',
        'weight' : 'normal',
        'size'   : 16,
        }
fontlabel = {'family' : 'serif',
        #'color'  : 'black',
        'weight' : 'normal',
        'size'   : 16,
        }


from matplotlib.ticker import FormatStrFormatter
plt.rc('font', **font)

```

# Data

## Probe locations

### Reflective blanket + mulch + concrete curing blanket

<table style="width:100%">
  <tr>
    <th>HOBO</th>
    <th>Channel</th> 
    <th>Location</th>
    <th>Note</th>
  </tr>
  <tr>
    <td>934</td>
    <td>2</td>
    <td>46 cm</td>
    <td> </td>
  </tr>
  <tr>
    <td>811</td>
    <td>1</td>
    <td>30 cm</td>
    <td> below reflective</td>
  </tr>
  <tr>
    <td>934</td>
    <td>1</td>
    <td>10 cm</td>
    <td> </td>
  </tr>
  <tr>
    <td>811</td>
    <td>1</td>
    <td>0 cm</td>
    <td> on snow</td>
  </tr>
</table>







```python
import pandas as pd
from pandas import Series

sensor811 = pd.read_excel("../Insulation/SN 20103811 2018-07-16 13_18_09 -0400.xlsx",sheet_name="DATA",
                          header=1,date_parser=[1],names=["#","timestamp","probe1","probe2"])
sensor811 = sensor811.set_index('timestamp')
sensor811 = sensor811.drop("#",axis=1)
sensor934 = pd.read_excel("../Insulation/SN 20366934 2018-07-16 13_19_37 -0400.xlsx",sheet_name="DATA",
                          header=1,date_parser=[0],names=["#","timestamp","probe1","probe2"])
sensor934 = sensor934.set_index('timestamp')
sensor934 = sensor934.drop("#",axis=1)
# print(sensor811.timestamp)
ax = sensor811.plot()
sensor934.plot(ax = ax)

alt0 = np.copy(sensor811.probe2)
alt30 = np.copy(sensor811.probe1)
alt10 = np.copy(sensor934.probe1)
alt46 = np.copy(sensor934.probe2)
```


### Reflective blanket + open cell foam + concrete curing blanket

<table style="width:100%">
  <tr>
    <th>HOBO</th>
    <th>Channel</th> 
    <th>Location</th>
    <th>Note</th>
  </tr>
  <tr>
    <td>930</td>
    <td>1</td>
    <td>46 cm</td>
    <td> </td>
  </tr>
  <tr>
    <td>930</td>
    <td>2</td>
    <td>30 cm</td>
    <td> below reflective</td>
  </tr>
  <tr>
    <td>812</td>
    <td>2</td>
    <td>10 cm</td>
    <td> </td>
  </tr>
  <tr>
    <td>812</td>
    <td>1</td>
    <td>0 cm</td>
    <td> on snow</td>
  </tr>
</table>



```python
import pandas as pd
from pandas import Series

sensor929 = pd.read_excel("../Insulation/SN 20366929 2018-07-16 13_18_44 -0400.xlsx",sheet_name="DATA",
                          header=1,date_parser=[1],names=["#","timestamp","probe1","probe2"])
sensor929 = sensor929.set_index('timestamp')
sensor929 = sensor929.drop("#",axis=1)
sensor930 = pd.read_excel("../Insulation/SN 20366930 2018-07-16 13_19_12 -0400.xlsx",sheet_name="DATA",
                          header=1,date_parser=[0],names=["#","timestamp","probe1","probe2"])
sensor930 = sensor930.set_index('timestamp')
sensor930 = sensor930.drop("#",axis=1)
# print(sensor811.timestamp)
ax = sensor929.plot()
sensor930.plot(ax = ax)

alt0 = np.copy(sensor929.probe2)
alt10 = np.copy(sensor929.probe1)
alt30 = np.copy(sensor930.probe1)
alt46 = np.copy(sensor930.probe2)
```


```python
tstart = sensor811.index[0]
print(tstart)
n = sensor811.shape[0]
print(n)
tend = sensor811.index[n-1]
print(tend)
print(pd.Timedelta(tend-tstart))

```

    2018-06-24 11:53:58
    6353
    2018-07-16 13:13:58
    22 days 01:20:00


# Time advancement




```python
dt = 5*60.

dt_2 = dt / 2
oneday = 24*3600
Sim_time = 22*oneday

Nt = int(Sim_time/dt)
print(oneday)
```

    86400


# Grid


```python
def create_grid(nz,lz,dzmin):
    z,gz = numtools.stretched_mesh(lz,nz,dzmin,28)
    
    dz = numtools.nufd(z)
    a_metrics = np.zeros(nz-2)
    c_metrics = np.zeros(nz-2)
    a_metrics[:] = 1. / ((z[1:-1] - z[0:-2])*(z[2:]-z[0:-2]))
    c_metrics[:] = 1. / ((z[2:] - z[1:-1])*(z[2:]-z[0:-2]))
    return z,dz,a_metrics,c_metrics
nz = 33
lz = 0.3
dzmin = lz/1000
z,dz,a_metrics,c_metrics = create_grid(nz,lz,dzmin)
plt.plot(z)
plt.show()
Dz = numtools.nufd(z)

```

# Simulation


```python
from scipy.interpolate import interp1d
z10 = np.array([0.1])

```


```python
def temperature_initialization(z,T_bot,T_top):
    return T_bot + (T_top - T_bot)*(z/z[-1]) 
    
```

## Constant Diffusivity Code


```python
nz = 65
lz = 0.2
dzmin = lz/1000
z,dz,a_metrics,c_metrics = create_grid(nz,lz,dzmin)
print(z[0],z[-1],lz)
start_time = 0*oneday
end_time = 22*oneday
alpha = 1e-5
a_rhs = np.zeros(nz-2)
b_rhs = np.zeros(nz-2)
c_rhs = np.zeros(nz-2)
a_lhs = np.zeros(nz-2)
b_lhs = np.zeros(nz-2)
c_lhs = np.zeros(nz-2)


T10 = np.zeros(Nt+1)


def bottom_T(i):
    global alt0
    return alt0[i]
def top_T(i):
    global alt30
    return alt30[i]
alpha_total_all = alpha*np.ones(nz)
a_lhs, b_lhs, c_lhs = numtools.diffusion_matrix_coefficients(alpha_total_all,
                                                                 dt_2,a_metrics,c_metrics)
istart = int(start_time/dt)
i = istart
T = temperature_initialization(z,alt0[i],alt30[i])
isample = 0
Sim_t = start_time
#     T[0] = bottom_T(i)
#     T[-1] = top_T(i)
if end_time > Sim_time:
    stop
while Sim_t < end_time:
#     print(Sim_t)
    Sim_t += dt
    a_rhs = np.copy(a_lhs)
    b_rhs = np.copy(b_lhs)
    c_rhs = np.copy(c_lhs)

    T_old = np.copy(T)
    T[0] = bottom_T(i)
    T[-1] = top_T(i)

    error = np.inf
    Titold = np.copy(T)
    j = 0
#     while error > 1e-5 and j < 100:
    alpha_total_all = alpha*np.ones(nz)
    a_lhs, b_lhs, c_lhs = numtools.diffusion_matrix_coefficients(alpha_total_all,
                                                                 dt_2,a_metrics,c_metrics)
    d = numtools.rhs_T(a_rhs, b_rhs, c_rhs, a_lhs, c_lhs, T_old, T,dt_2,dt)
    T[1:-1] = numtools.TDMAsolver(-a_lhs, 1-b_lhs , -c_lhs, d)
#         error = np.linalg.norm(T-Titold,2)/np.mean(T)
#         Titold = np.copy(T)
#         j += 1
#             print(j,error)
    i += 1
    f = interp1d(z,T)
    T10[i] = f(z10)

plt.plot(T10)
plt.plot(alt10)
plt.show()
plt.plot(z,T)
plt.show()
```


```python

def cst_diffusivity_sim(alpha,L_ccb,start_time,end_time):
    global z,dt, dt_2, Sim_time, Nt, z10, a_metrics, c_metrics,alt0,alt30,Dz

    a_rhs = np.zeros(nz-2)
    b_rhs = np.zeros(nz-2)
    c_rhs = np.zeros(nz-2)
    a_lhs = np.zeros(nz-2)
    b_lhs = np.zeros(nz-2)
    c_lhs = np.zeros(nz-2)

    Tprofiles = np.zeros((nz,Nt+1))
    T10 = np.zeros(Nt+1)


    def bottom_T(i):
        global alt0
        return alt0[i]
    def top_T(i):
        global alt30
        return alt30[i]
    mask_ccb = np.where(z < L_ccb)
    alpha_total_all = alpha[1]*np.ones(nz)
    alpha_total_all[mask_ccb] = alpha[0]
    istart = int(start_time/dt)
    i = istart
    T = temperature_initialization(z,alt0[i],alt30[i])
    a_lhs, b_lhs, c_lhs = numtools.diffusion_matrix_coefficients(alpha_total_all,dt_2,a_metrics,c_metrics)
    isample = 0
    Sim_t = start_time
#     T[0] = bottom_T(i)
#     T[-1] = top_T(i)
    if end_time > Sim_time:
        stop
    while Sim_t < end_time:

        Sim_t += dt
        a_rhs = np.copy(a_lhs)
        b_rhs = np.copy(b_lhs)
        c_rhs = np.copy(c_lhs)

        T_old = np.copy(T)
        T[0] = bottom_T(i)
        T[-1] = top_T(i)

        error = np.inf
        Titold = np.copy(T)
        j = 0
        while error > 1e-5 and j < 100:
            a_lhs, b_lhs, c_lhs = numtools.diffusion_matrix_coefficients(alpha_total_all,dt_2,a_metrics,c_metrics)
            d = numtools.rhs_T(a_rhs, b_rhs, c_rhs, a_lhs, c_lhs, T_old, T,dt_2,dt)
            T[1:-1] = numtools.TDMAsolver(-a_lhs, 1-b_lhs , -c_lhs, d)
            error = np.linalg.norm(T-Titold,2)/np.mean(T)
            Titold = np.copy(T)
            j += 1
#             print(j,error)
        i += 1
        f = interp1d(z,T)
        T10[i] = f(z10)
        Tprofiles[:,i] = np.copy(T)
#         print(i,Sim_t/Sim_time*100)
#         clear_output(wait=True)
    return T10,Tprofiles

```


```python
alpha_ccb = 1e-7
alpha_mulch = 3e-7
L_ccb = 0.0254
x = np.array([alpha_ccb, alpha_mulch])
nz = 33
lz = 0.2
dzmin = lz/100
start_time = 0*oneday
end_time = 22*oneday
istart = int(start_time/dt)
iend = int(end_time/dt)
z,dz,a_metrics,c_metrics = create_grid(nz,lz,dzmin)
T10,Tprofiles = cst_diffusivity_sim(x,L_ccb,start_time,end_time)
plt.plot(T10)
plt.plot(alt10)
plt.show()
```


```python
nzarray = np.array([17,33,65,129])
for nz in nzarray:
    dzmin = lz/100
    z,dz,a_metrics,c_metrics = create_grid(nz,lz,dzmin)
    T10,Tprofiles = cst_diffusivity_sim(x,L_ccb,start_time,end_time)
    plt.plot(T10)
plt.show()
```


```python
nz = 33
lz = 0.2
dzmin = lz/1000
z,dz,a_metrics,c_metrics = create_grid(nz,lz,dzmin)
```


```python
n = 5
alpha_mulch = np.logspace(-7,-6,n)
alpha_ccb = np.logspace(-7.5,-6.5,n)
X,Y = np.meshgrid(alpha_ccb,alpha_mulch)
error_alpha_cst = np.zeros((n,n))
i = 0
# T10_all = np.zeros((Nt+1,n))
i = 0
istart = int(start_time/dt)
iend =int(end_time/dt)
for alpha_m in alpha_mulch:
    j = 0
    for alpha_c in alpha_ccb:
        x = np.array([alpha_c,alpha_m])
        T10,Tprofiles = cst_diffusivity_sim(x, L_ccb,start_time,end_time)
#     T10_all[:,i] = np.copy(T10)
        
        error_alpha_cst[i,j] = np.std(T10[istart+1:iend] - alt10[istart+1:iend])
        print(i,j)
        j += 1
        clear_output(wait=True)
    i += 1
plt.contourf(X,Y,error_alpha_cst,100)
plt.xscale('log')
plt.yscale('log')
plt.show()
```


```python
plt.contourf(X,Y,error_alpha_cst,100,cmap = 'jet')
plt.xscale('log')
plt.yscale('log')
plt.show()

```


```python
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cbook
from matplotlib import cm
from matplotlib.colors import LightSource

fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))

ls = LightSource(270, 45)
# To use a custom hillshading mode, override the built-in shading and pass
# in the rgb colors of the shaded surface calculated from "shade".
rgb = ls.shade(error_alpha_cst, cmap=cm.gist_earth, vert_exag=0.1, blend_mode='soft')
surf = ax.plot_surface(X,Y,error_alpha_cst, rstride=1, cstride=1, facecolors=rgb,
                       linewidth=0, antialiased=False, shade=False)
ax.view_init(45, 60)
plt.show()
```

$$
J(\alpha)=\min_\alpha\left(\left\Vert T_{sim}(z=10cm,\alpha)-T_{HOBO}(z=10cm)\right\Vert\right)
$$



```python
start_time = 10*oneday
end_time = 20*oneday
istart = int(start_time/dt)
iend = int(end_time/dt)
def cost_function(x):
    global z,dt, dt_2, Sim_time, Nt, z10, a_metrics, c_metrics,alt0 \
    ,alt30,alt10,istart,iend,start_time,end_time,L_ccb
    T10,Tprofiles = cst_diffusivity_sim(x,L_ccb,start_time,end_time)
    return np.std(T10[istart+1:iend] - alt10[istart+1:iend])
from scipy import optimize
x0 = [5e-7,3e-7]
x_opt = optimize.minimize(cost_function,x0, method='Nelder-Mead',
#                           jac='2-point',
#                           bounds=((1e-8,1e-6),(1e-7,1e-6)), 
                     options={'maxiter':50,'disp':True})
print(x_opt.nit,x_opt.message,x_opt.x)
```

    /Users/dubief/gitrepo/snow-farming/Modeling/numericaltools.py:42: RuntimeWarning: overflow encountered in double_scalars
      xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]
    /Users/dubief/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:53: RuntimeWarning: invalid value encountered in double_scalars
    /Users/dubief/gitrepo/snow-farming/Modeling/numericaltools.py:82: RuntimeWarning: overflow encountered in multiply
      d = a_rhs*T_old_all[:-2] + (b_rhs + 1)*T_old_all[1:-1] + c_rhs*T_old_all[2:]
    /Users/dubief/gitrepo/snow-farming/Modeling/numericaltools.py:82: RuntimeWarning: invalid value encountered in add
      d = a_rhs*T_old_all[:-2] + (b_rhs + 1)*T_old_all[1:-1] + c_rhs*T_old_all[2:]


    Optimization terminated successfully.
             Current function value: 0.286843
             Iterations: 24
             Function evaluations: 45
    24 Optimization terminated successfully. [1.30112648e-07 3.85262403e-07]



```python
# x = np.array([1e-7, 5e-7])
x = x_opt.x
print(x)
T10,Tprofiles = cst_diffusivity_sim(x,L_ccb,0,Sim_time)
plt.plot(T10[istart+1:iend])
plt.plot(alt10[istart+1:iend])
plt.show()
plt.plot(T10)
plt.plot(alt10)
plt.show()
diff = (T10[int(oneday/dt):iend]-alt10[int(oneday/dt):iend])
DT = (alt30[int(oneday/dt)]-alt0[int(oneday/dt)])
norm = (T10[int(oneday/dt):iend]-alt10[int(oneday/dt):iend])/(alt30[int(oneday/dt)]-alt0[int(oneday/dt)])
plt.plot(norm)
plt.show()
```


```python
for i in range(0,6000,100):
    plt.plot(z,Tprofiles[:,i])
plt.show()
```


```python
import pandas as pd
from pandas import Series

weather1 = pd.read_csv("../Excel/Weather_KVTCRAFT2_2017-06-11_2018-03-20.csv", delimiter=",",header=0,date_parser=[1])
weather2 = pd.read_csv("../Excel/Weather_KVTCRAFT2_2018-03-23_2018-10-16.csv", delimiter=",",header=0,date_parser=[1])
weather = [weather1, weather2]
weather = pd.concat(weather)
weather['time'] = pd.to_timedelta(weather['time'].astype(str))
weather['date'] = pd.to_datetime(weather['date'])
weather['date'] = weather['date'] + weather['time']
weather['date'] = pd.to_datetime(weather['date'])

# weather.set_index('date')
weather.columns
```




    Index(['date', 'time', 'temperature', 'dewpoint', 'humidity', 'wind_speed',
           'wind_gust_speed', 'wind_dir_degrees', 'wind_dir', 'pressure',
           'windchill', 'heatindex', 'precip', 'precip_rate', 'precip_1hr',
           'precip_today', 'solarradiation', 'uv_index', 'temperature_indoor',
           'humidity_indoor', 'software_type', 'event', 'conditions'],
          dtype='object')




```python
ax = weather.plot(x='date', y = 'solarradiation')
ax.set_xlim(sensor811.index[0],sensor811.index[-1])
```


```python
elaspsedtime_mess = np.array(weather['date'].astype('str'))
elaspedtime=np.array([])
irradiation = np.array(weather['solarradiation'])
from datetime import datetime
for i in range(2): #range(len(elaspsedtime_mess)):
    strdate = elaspsedtime_mess[i]
    year = strdate[:4]
    month = strdate[5:7]
    day = strdate[8:11]
    
    hour = strdate[10:13]
    minute = strdate[14:16]
    second = strdate[15:16]
    print(year,month,day,hour,minute,second)
print(elpasedtime[1])
print(irradiation[:5])
```

    2017 06 11   00 01 1
    2017 06 11   00 31 1
    2017-06-11T00:31:00.000000000
    [0. 0. 0. 0. 0.]



```python
plt.plot(dzT_top)
plt.plot(dzT_bottom*0.08)
plt.xlim(0,6000)
plt.ylim(-10,100)
plt.grid(True)
```

# Volume Loss

Equation of melt, where $V$ is the snow volume, $rho_{snow}$ and $h_{L,snow}$ are the snow density and latent heat, respectively.

$$
\rho_{snow}h_{L,snow}\frac{dV}{dt}=
\begin{cases}
-q_{outside} & \text{when }q_{outside}>0\\
0 &\text{otherwise}
\end{cases}
$$


Simplest model for volume loss:



$$
\frac{dV}{dt}=-\frac{q}{\rho_{snow}h_{L,snow}} = -\frac{1}{\rho_{snow}h_{L,snow}}\left(\frac{T_{atm}-T_{snow}}{R}+q_{sky}+q_{sun}\right)
$$


```python
Volumes = pd.read_excel("../Excel/Volumes_snow.xlsx",
                    header = 0,
                    names = ["Date","Site 1","Site 2"])

# Volumes = Volumes.dropna()

Volumes.set_index('Date')
Volumes['Site 1']
Volumes1 = Volumes[1:]
```


```python
Volumes1.plot(x='Date',y='Site 2')
```


```python
Volumes1['date_ordinal'] = pd.to_datetime(Volumes1['Date']).apply(lambda date: date.toordinal())
Volumes1
```

    /Users/dubief/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
      """Entry point for launching an IPython kernel.





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Date</th>
      <th>Site 1</th>
      <th>Site 2</th>
      <th>date_ordinal</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>1</th>
      <td>2018-04-21</td>
      <td>249.07</td>
      <td>251.17</td>
      <td>736805</td>
    </tr>
    <tr>
      <th>2</th>
      <td>2018-05-12</td>
      <td>216.50</td>
      <td>226.13</td>
      <td>736826</td>
    </tr>
    <tr>
      <th>3</th>
      <td>2018-05-28</td>
      <td>195.50</td>
      <td>208.10</td>
      <td>736842</td>
    </tr>
    <tr>
      <th>4</th>
      <td>2018-06-10</td>
      <td>177.00</td>
      <td>187.70</td>
      <td>736855</td>
    </tr>
    <tr>
      <th>5</th>
      <td>2018-06-20</td>
      <td>157.23</td>
      <td>170.13</td>
      <td>736865</td>
    </tr>
    <tr>
      <th>6</th>
      <td>2018-06-30</td>
      <td>139.10</td>
      <td>150.90</td>
      <td>736875</td>
    </tr>
    <tr>
      <th>7</th>
      <td>2018-07-03</td>
      <td>83.03</td>
      <td>142.47</td>
      <td>736878</td>
    </tr>
    <tr>
      <th>8</th>
      <td>2018-07-10</td>
      <td>70.07</td>
      <td>130.33</td>
      <td>736885</td>
    </tr>
    <tr>
      <th>9</th>
      <td>2018-07-21</td>
      <td>61.67</td>
      <td>116.10</td>
      <td>736896</td>
    </tr>
    <tr>
      <th>10</th>
      <td>2018-07-31</td>
      <td>50.27</td>
      <td>100.83</td>
      <td>736906</td>
    </tr>
    <tr>
      <th>11</th>
      <td>2018-08-09</td>
      <td>42.83</td>
      <td>91.20</td>
      <td>736915</td>
    </tr>
    <tr>
      <th>12</th>
      <td>2018-08-23</td>
      <td>33.73</td>
      <td>75.83</td>
      <td>736929</td>
    </tr>
    <tr>
      <th>13</th>
      <td>2018-09-09</td>
      <td>28.77</td>
      <td>60.87</td>
      <td>736946</td>
    </tr>
  </tbody>
</table>
</div>




```python
ax = sns.regplot(
    data=Volumes1,
    x='date_ordinal',
    y='Site 2',
)
# Tighten up the axes for prettiness
ax.set_xlim(Volumes1['date_ordinal'].min() - 1, Volumes1['date_ordinal'].max() + 1)
ax.set_ylim(0, Volumes1['Site 2'].max() + 1)
```


```python
ax.set_xlabel('Date')
new_labels = [Date.fromordinal(int(item)) for item in ax.get_xticks()]
ax.set_xticklabels(new_labels)
```


```python
import statsmodels.api as sm
y = Volumes1["Site 2"][:]
X = Volumes1["date_ordinal"]
model = sm.OLS(y,X)
results = model.fit()
print(results.summary())
```

                                OLS Regression Results                            
    ==============================================================================
    Dep. Variable:                 Site 2   R-squared:                       0.870
    Model:                            OLS   Adj. R-squared:                  0.859
    Method:                 Least Squares   F-statistic:                     80.23
    Date:                Wed, 10 Apr 2019   Prob (F-statistic):           1.16e-06
    Time:                        10:37:29   Log-Likelihood:                -70.976
    No. Observations:                  13   AIC:                             144.0
    Df Residuals:                      12   BIC:                             144.5
    Df Model:                           1                                         
    Covariance Type:            nonrobust                                         
    ================================================================================
                       coef    std err          t      P>|t|      [0.025      0.975]
    --------------------------------------------------------------------------------
    date_ordinal     0.0002   2.23e-05      8.957      0.000       0.000       0.000
    ==============================================================================
    Omnibus:                        0.846   Durbin-Watson:                   0.077
    Prob(Omnibus):                  0.655   Jarque-Bera (JB):                0.705
    Skew:                           0.260   Prob(JB):                        0.703
    Kurtosis:                       1.984   Cond. No.                         1.00
    ==============================================================================
    
    Warnings:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.


    /Users/dubief/anaconda3/lib/python3.6/site-packages/scipy/stats/stats.py:1416: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=13
      "anyway, n=%i" % int(n))



```python
print('Parameters: ', results.params)
print('R2: ', results.rsquared)

Volumes1["LinReg"] = 0.0002
```

    Parameters:  date_ordinal    0.0002
    dtype: float64
    R2:  0.8698962333612299



```python

def cst_source_diffusivity_sim(x,L_ccb,start_time,end_time):
    global z,dt, dt_2, Sim_time, Nt, z10, a_metrics, c_metrics,alt0,alt30

    a_rhs = np.zeros(nz-2)
    b_rhs = np.zeros(nz-2)
    c_rhs = np.zeros(nz-2)
    a_lhs = np.zeros(nz-2)
    b_lhs = np.zeros(nz-2)
    c_lhs = np.zeros(nz-2)

    
    T10 = np.zeros(Nt+1)


    def bottom_T(i):
        global alt0
        return alt0[i]
    def top_T(i):
        global alt30
        return alt30[i]
    mask_ccb = np.where(z < L_ccb)
    alpha_total_all = x[1]*np.ones(nz)
    alpha_total_all[mask_ccb] = x[0]
    istart = int(start_time/dt)
    i = istart
    T = temperature_initialization(z,alt0[i],alt30[i])
    a_lhs, b_lhs, c_lhs = numtools.diffusion_matrix_coefficients(alpha_total_all,dt_2,a_metrics,c_metrics)
    isample = 0
    Sim_t = start_time
#     T[0] = bottom_T(i)
#     T[-1] = top_T(i)
    if end_time > Sim_time:
        stop
    while Sim_t < end_time:

        Sim_t += dt
        a_rhs = np.copy(a_lhs)
        b_rhs = np.copy(b_lhs)
        c_rhs = np.copy(c_lhs)

        T_old = np.copy(T)
        T[0] = bottom_T(i)
        T[-1] = top_T(i)
        DeltaT = T[-1] - T[0]
        error = np.inf
        Titold = np.copy(T)
        j = 0
        while error > 1e-5 and j < 100:
            a_lhs, b_lhs, c_lhs = numtools.diffusion_matrix_coefficients(alpha_total_all,dt_2,a_metrics,c_metrics)
            d = numtools.rhs_T(a_rhs, b_rhs, c_rhs, a_lhs, c_lhs, T_old, T,dt_2,dt)
            d += x[2]*DeltaT
            T[1:-1] = numtools.TDMAsolver(-a_lhs, 1-b_lhs , -c_lhs, d)
            error = np.linalg.norm(T-Titold,2)/np.mean(T)
            Titold = np.copy(T)
            j += 1
#             print(j,error)
        i += 1
        f = interp1d(z,T)
        T10[i] = f(z10)
#         print(i,Sim_t/Sim_time*100)
#         clear_output(wait=True)
    return T10
```


```python
x0 = [1.58659189e-07, 3.53243639e-07, 7.06084974e-04]
T10 = cst_source_diffusivity_sim(x0,L_ccb, 0,Sim_time)
plt.plot(T10[istart+1:iend])
plt.plot(alt10[istart+1:iend])
plt.show()
plt.plot(T10)
plt.plot(alt10)
plt.show()
```


```python
diff = (T10[int(oneday/dt):iend]-alt10[int(oneday/dt):iend])
DT = (alt30[int(oneday/dt):iend]-alt0[int(oneday/dt):iend])
T_f = (alt30[int(oneday/dt):iend]+alt0[int(oneday/dt):iend])/2
norm = diff / DT
plt.plot(norm)
plt.show()
plt.plot(DT,diff)
plt.show()
plt.plot(T_f,diff)
plt.show()
```


```python
def cost_function_cst_source(x):
    global z,dt, dt_2, Sim_time, Nt, z10, a_metrics, c_metrics,alt0 \
    ,alt30,alt10,istart,iend,start_time,end_time
    T10 = cst_source_diffusivity_sim(x,L_ccb,start_time,end_time)
    return np.std(T10[istart+1:iend] - alt10[istart+1:iend])
from scipy import optimize
x0 = [1.58659189e-07, 3.53243639e-07, 7.06084974e-04]
x_opt = optimize.minimize(cost_function_cst_source,x0, method='Nelder-Mead',
                     options={'maxiter':200,'disp':True})
print(x_opt.nit,x_opt.message,x_opt.x)
```


```python

T10 = cst_source_diffusivity_sim(x_opt.x,L_ccb,0,Sim_time)
plt.plot(T10[istart+1:iend])
plt.plot(alt10[istart+1:iend])
plt.show()
plt.plot(T10)
plt.plot(alt10)
plt.show()
```


```python
def lin_diffusivity_sim(x,start_time,end_time):
    global z,dt, dt_2, Sim_time, Nt, z10, a_metrics, c_metrics,alt0,alt30
    
    def alpha_fct(x,T):
        theta = (T-5)
        mask = np.where(theta < 0)
        theta[mask] = 0
        
        return x[0]*(1+x[1]*theta)

    a_rhs = np.zeros(nz-2)
    b_rhs = np.zeros(nz-2)
    c_rhs = np.zeros(nz-2)
    a_lhs = np.zeros(nz-2)
    b_lhs = np.zeros(nz-2)
    c_lhs = np.zeros(nz-2)

    T = temperature_initialization(z,alt0[0],alt30[0])
    T10 = np.zeros(Nt+1)


    def bottom_T(i):
        global alt0
        return alt0[i]
    def top_T(i):
        global alt30
        return alt30[i]
    istart = int(start_time/dt)
    i = istart
    T = temperature_initialization(z,alt0[i],alt30[i])
    isample = 0
    Sim_t = start_time
#     T[0] = bottom_T(i)
#     T[-1] = top_T(i)
    if end_time > Sim_time:
        stop
    while Sim_t < end_time:

        Sim_t += dt
        a_rhs = np.copy(a_lhs)
        b_rhs = np.copy(b_lhs)
        c_rhs = np.copy(c_lhs)

        T_old = np.copy(T)
        T[0] = bottom_T(i)
        T[-1] = top_T(i)

        error = np.inf
        Titold = np.copy(T)
        j = 0
        while error > 1e-5 and j < 100:
            alpha_total_all = alpha_fct(x,T)
            a_lhs, b_lhs, c_lhs = numtools.diffusion_matrix_coefficients(alpha_total_all,dt_2,a_metrics,c_metrics)
            d = numtools.rhs_T(a_rhs, b_rhs, c_rhs, a_lhs, c_lhs, T_old, T,dt_2,dt)
            T[1:-1] = numtools.TDMAsolver(-a_lhs, 1-b_lhs , -c_lhs, d)
            error = np.linalg.norm(T-Titold,2)/np.mean(T)
            Titold = np.copy(T)
            j += 1
#             print(j,error)
        i += 1
        f = interp1d(z,T)
        T10[i] = f(z10)
#         print(i,Sim_t/Sim_time*100)
#         clear_output(wait=True)
    return T10


```


```python
x = [2e-7, 1]
T10 = lin_diffusivity_sim(x,start_time,end_time)
plt.plot(T10)
plt.plot(alt10)
plt.show()
```


```python
start_time = 10*oneday
end_time = 12*oneday
istart = int(start_time/dt)
iend = int(end_time/dt)
from scipy import optimize
def cost_function_lin(x):
    global z,dt, dt_2, Sim_time, Nt, z10, a_metrics, c_metrics,alt0,alt30,alt10,istart,iend,start_time,end_time 
    T10 = lin_diffusivity_sim(x,start_time,end_time)
    return np.std(T10[istart+1:iend] - alt10[istart+1:iend])
x0 = [1.3e-6, -0.04]
x_opt = optimize.minimize(cost_function_lin,x0, method='Nelder-Mead',
                     options={'maxiter':200,'disp':True})
print(x_opt.nit,x_opt.message,x_opt.x)
```


```python

```


```python
T10 = lin_diffusivity_sim(x_opt.x,oneday,Sim_time)
plt.plot(T10)
plt.plot(alt10)
plt.show()
plt.plot(T10[2*int(oneday/dt):iend]-alt10[2*int(oneday/dt):iend])
plt.show()
x0_lin = np.copy(x_opt.x)
```


```python

```


```python
def lin_S_diffusivity_sim(x,start_time,end_time):
    global z,dt, dt_2, Sim_time, Nt, z10, a_metrics, c_metrics,alt0,alt30
    
    def alpha_fct(x,T):
        theta = (T-5)
        mask = np.where(theta < 0)
        theta[mask] = 0
        
        return x[0]*(1+x[1]*theta)

    a_rhs = np.zeros(nz-2)
    b_rhs = np.zeros(nz-2)
    c_rhs = np.zeros(nz-2)
    a_lhs = np.zeros(nz-2)
    b_lhs = np.zeros(nz-2)
    c_lhs = np.zeros(nz-2)

    T = temperature_initialization(z,alt0[0],alt30[0])
    T10 = np.zeros(Nt+1)


    def bottom_T(i):
        global alt0
        return alt0[i]
    def top_T(i):
        global alt30
        return alt30[i]
    istart = int(start_time/dt)
    i = istart
    T = temperature_initialization(z,alt0[i],alt30[i])
    isample = 0
    Sim_t = start_time
#     T[0] = bottom_T(i)
#     T[-1] = top_T(i)
    if end_time > Sim_time:
        stop
    while Sim_t < end_time:

        Sim_t += dt
        a_rhs = np.copy(a_lhs)
        b_rhs = np.copy(b_lhs)
        c_rhs = np.copy(c_lhs)

        T_old = np.copy(T)
        T[0] = bottom_T(i)
        T[-1] = top_T(i)

        error = np.inf
        Titold = np.copy(T)
        j = 0
        while error > 1e-5 and j < 100:
            alpha_total_all = alpha_fct(x,T)
            a_lhs, b_lhs, c_lhs = numtools.diffusion_matrix_coefficients(alpha_total_all,dt_2,a_metrics,c_metrics)
            d = numtools.rhs_T(a_rhs, b_rhs, c_rhs, a_lhs, c_lhs, T_old, T,dt_2,dt)
            T[1:-1] = numtools.TDMAsolver(-a_lhs, 1-b_lhs , -c_lhs, d)
            error = np.linalg.norm(T-Titold,2)/np.mean(T)
            Titold = np.copy(T)
            j += 1
#             print(j,error)
        i += 1
        f = interp1d(z,T)
        T10[i] = f(z10)
#         print(i,Sim_t/Sim_time*100)
#         clear_output(wait=True)
    return T10


```


```python

```

c406c13d62e1e3dfbfa548d19527e092127e2852
ipynb
Jupyter Notebook
A01-Grid-Generation/.ipynb_checkpoints/grid-generation-checkpoint.ipynb
yvesdubief/numericalmethods
725c3f08159e3ce62f49035c74a3725401cce9df
MIT
2020-10-18T19:40:22.000Z
2020-10-18T19:40:22.000Z
A01-Grid-Generation/.ipynb_checkpoints/grid-generation-checkpoint.ipynb
fuqianyin/numericalmethods
e8ac8e4c378fa1c6ed9cfd67a89bcc12567f658d
MIT
A01-Grid-Generation/.ipynb_checkpoints/grid-generation-checkpoint.ipynb
fuqianyin/numericalmethods
e8ac8e4c378fa1c6ed9cfd67a89bcc12567f658d
MIT
2019-10-11T15:24:13.000Z
2019-10-11T15:24:13.000Z
Qwen/Qwen-72B
1. YES
2. YES