# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ### This notebook compares Turn Point Runs using 1000, 2000 and 3000 particles

# +
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import netCDF4 as nc
import seaborn as sns
import matplotlib.colors as mcolors
from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.patches as mpatches
import cmocean as cm
import glob
import os
import xarray as xr
import datetime
from salishsea_tools import viz_tools, tidetools, geo_tools, gsw_calls, wind_tools
import pickle
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

from sys import platform

# %matplotlib inline

if platform == 'darwin':
    bathy_dir = '/Users/rmueller/Data/SalishSeaCast/grid/'
    results_dir = '/Users/rmueller/Projects/MIDOSS/results/nbr_partic_test/'
    out_dir = '/Users/rmueller/Documents/Presentations/Graphics/MIDOSS/SOILED/nbr_partic/'
else:
    display('Update file paths for oceans machine')
# -

# define number of particles and the compute time for each run
nbr_partic = [1000,2000,3000]
compute_time = [102, 104, 102]
cpu_time = [5934.55, 6025.39, 5670.18]

tp_akns={} # this is the way to define a new dictionary
thick2D={} # one dictionary is never enough
thick2D_t={} # every house ought to have at least three!
for nparticles in nbr_partic:
    tp_akns[nparticles]  = xr.open_dataset(f'{results_dir}Lagrangian_AKNS_crude_TP_01aug17-08aug17_akns{nparticles}.nc')
    thick2D[nparticles]  = tp_akns[nparticles].Thickness_2D
    thick2D_t[nparticles] = thick2D[nparticles].sum(dim='time')

# ### Use sum of thickness 2D to help discern whether the mass discrepency is an output error or a mass calculation error

[nt,ny,nx] = thick2D[1000].shape
np.sum(thick2D[1000][1,:,:].values)

# +
thick2D_timeseries = {}

for nparticles in nbr_partic:
    
    thick2D_timeseries[nparticles] = np.array([])
   
    for time_index in range(nt):
        
        timeseries_value = np.sum(thick2D[nparticles][time_index,:,:].values)   
        # The first value of sum is negative so I explicitely calculate it here to make sure the results isn't a funky issue with the append function, which it isn't
        if time_index == 0:
            thick2D_timeseries[nparticles] = timeseries_value
        else:
            thick2D_timeseries[nparticles] = np.append(thick2D_timeseries[nparticles], timeseries_value)
# -



# +
fs = 20
fig = plt.figure(figsize=(20,20))
ax1 = fig.add_subplot(131)
ax2 = fig.add_subplot(132)
ax3 = fig.add_subplot(133)
# convert xarray into numpy using ".values" in order to gain access to different visualization tools
mappable = ax1.pcolormesh(thick2D_t[1000].values, vmin = 0, vmax = 40, cmap = cm.cm.balance)
mappable = ax2.pcolormesh(thick2D_t[2000].values, vmin = 0, vmax = 40, cmap = cm.cm.balance)
mappable = ax3.pcolormesh(thick2D_t[3000].values, vmin = 0, vmax = 40, cmap = cm.cm.balance)

# add land mask to ax1 and ax2
viz_tools.plot_land_mask(ax1,'/Users/rmueller/Projects/MIDOSS/MIDOSS-MOHID-grid/AfterNEMOBathy201702.nc', color = 'burlywood')
viz_tools.plot_land_mask(ax2,'/Users/rmueller/Projects/MIDOSS/MIDOSS-MOHID-grid/AfterNEMOBathy201702.nc', color = 'burlywood')
viz_tools.plot_land_mask(ax3,'/Users/rmueller/Projects/MIDOSS/MIDOSS-MOHID-grid/AfterNEMOBathy201702.nc', color = 'burlywood')

ax1.set_title(r'1000 ($\sum$ thickness =' + f'{thick2D_timeseries[1000][nt-1]:5.1f})', fontsize = fs)
ax2.set_title(r'2000 ($\sum$ thickness =' + f'{thick2D_timeseries[2000][nt-1]:5.1f})', fontsize = fs)
ax3.set_title(r'3000 ($\sum$ thickness =' + f'{thick2D_timeseries[3000][nt-1]:5.1f})', fontsize = fs)

ax1.set_ylim(200,400)
ax2.set_ylim(200,400)
ax3.set_ylim(200,400)

ax1.set_xlim(50,350)
ax2.set_xlim(50,350)
ax3.set_xlim(50,350)

ax1.set_aspect(aspect=1)
ax2.set_aspect(aspect=1)
ax3.set_aspect(aspect=1)
# -

# ### Check to make sure that the mass of oil in these three cases is the same

# +
#### Load header information

with open(f'{results_dir}resOilOutput_1000.sro', 'r') as the_file:
    all_data = [line.strip() for line in the_file.readlines()]
    header = all_data[4]  
    
# Order header into list array by splitting up string
header_arr = []
header_arr = header.split(' ')

# Remove emtpy entries from list
header_arr = np.asarray([x for x in header_arr if x != ''])

# load mass balance area values
mass_bal = {}
mass_bal_infile = {}

for nparticles in nbr_partic:
    mass_bal_infile[nparticles] = f'{results_dir}resOilOutput_{nparticles}.sro'
    mass_bal[nparticles]        = np.genfromtxt(mass_bal_infile[nparticles], skip_header=6, skip_footer=4)

# -

header_arr

#number of time points and parameters
[nt,npoints] = mass_bal[1000].shape

# indices for mass balance (etc)
i_mevaporated = 15
i_mdispersed  = 18
i_mdissolved  = 24
i_mbio        = 37
i_voloilbeached = 8
i_volumebeached = 9
i_volumeoil     = 10
i_volume        = 11
i_vwatercontent = 34
i_msedimented   = 21
i_mwatercontent = 33
i_density   = 35
i_viscosity = 36 
i_area = 12
i_theorical_area = 13
# analyte mass (I thought this was dissolved but the numbers are more reflective of dispersed)
i_analytemass0 = 42
i_analytemass1 = 43
i_analytemass2 = 44
i_analytemass3 = 45
i_analytemass4 = 46
# biodegredation
i_bio0 = 47
i_bio1 = 48
i_bio2 = 49
i_bio3 = 50
i_bio4 = 51

# plot up volume 
plt.plot(mass_bal[1000][range(nt),i_volume],'b')
plt.plot(mass_bal[2000][range(nt),i_volume])
plt.plot(mass_bal[3000][range(nt),i_volume],'g')
plt.ylabel(header_arr[i_volume])
plt.xlabel('Hours after spill ')
plt.legend(['1000', '2000', '3000'])
plt.title('NBR_PARTIC comparison')

plot_mass = [i_mevaporated, i_mdispersed, i_mdissolved, i_mbio]
for imass in plot_mass: 
    for nparticles in nbr_partic:
        plt.plot(mass_bal[nparticles][range(nt),imass])
plt.ylabel('Mass (Tonnes)')
plt.xlabel('Time after oil release (hours) ')

plot_mass = [i_mevaporated, i_mdispersed, i_mdissolved, i_mbio]
for imass in plot_mass: 
    plt.figure()
    for nparticles in nbr_partic:  
        plt.plot(mass_bal[nparticles][range(nt),imass])
    plt.legend(['1000','2000','3000'])        
    plt.ylabel('Mass (Tonnes)')
    plt.xlabel('Time after oil release (hours) ')
    plt.title(header_arr[imass])

# ### Evaluate results from changing Mass used in MDispersedDT calculation from Me%Var%MassOil to Me%Var%MassINI

# +
nbr_partic = [1000,2000, 3000, 3001, 3002]

# load mass balance area values
mass_bal = {}
mass_bal_infile = {}

for nparticles in nbr_partic:
    if nparticles==3001:
        mass_bal_infile[nparticles] = f'{results_dir}resOilOutput_3000_fix1.sro'
        mass_bal[nparticles]        = np.genfromtxt(mass_bal_infile[nparticles], skip_header=6, skip_footer=4)
    elif nparticles==3002:
        mass_bal_infile[nparticles] = f'{results_dir}resOilOutput_3000_fix2.sro'
        mass_bal[nparticles]        = np.genfromtxt(mass_bal_infile[nparticles], skip_header=6, skip_footer=4)
    else:
        mass_bal_infile[nparticles] = f'{results_dir}resOilOutput_{nparticles}.sro'
        mass_bal[nparticles]        = np.genfromtxt(mass_bal_infile[nparticles], skip_header=6, skip_footer=4)
    


# -

plot_mass = [i_mevaporated, i_mdispersed, i_mdissolved, i_mbio]
for imass in plot_mass: 
    plt.figure()
    for nparticles in nbr_partic:  
        plt.plot(mass_bal[nparticles][range(nt),imass])
    plt.legend(['1000','2000','3000','3000(MassINI)','3000(bio mass fix)'])        
    plt.ylabel('Mass (Tonnes)')
    plt.xlabel('Time after oil release (hours) ')
    plt.title(header_arr[imass])

mass_bal_infile[nparticles]

# #### Test whether mass issue is an output or whether the 2D thickness values also show different mass values 

# +
## calculate the increase in surface mass when comparing NBR_PARTIC = 2000 and 3000 to 1000.  
# -

# The first value of this sum is negative, so I'm plotting subsequent values 
surfacemass_timeseries = {}
for nparticles in nbr_partic[0:3]:
    surfacemass_timeseries[nparticles] = 500 * 440 *  mass_bal[nparticles][range(nt),i_density] * thick2D_timeseries[nparticles] * 1e-9
    plt.plot(surfacemass_timeseries[nparticles][1:-1])
    display(f'first value < 0 for nbrpartic[{nparticles}]: {thick2D_timeseries[nparticles][0]} microns.  Plotting subsequent values.')
plt.ylabel('Integrated surface mass (tonnes)')
plt.xlabel('Time after oil release (hours) ')
plt.title('Spatially-integrated surface mass (tonnes) after first saved output (which is negative)')
plt.legend([f'1000 (t(end) = {surfacemass_timeseries[1000][nt-1]:3.1f})', 
          f'2000 (t(end) = {surfacemass_timeseries[2000][nt-1]:3.1f})',
          f'3000 (t(end) = {surfacemass_timeseries[3000][nt-1]:3.1f})'
           ])

# ### verify that the mass is comparable in the output from the version of NBR_PARTIC = 3000 pre- and post- bug fix to remove biodegredation mass (which is << than other mass losses; hence, the comparison ought to be comparable)

tp_akns[3002]  = xr.open_dataset(f'{results_dir}Lagrangian_AKNS_crude_TP_01aug17-08aug17_akns3000_fix2.nc')
thick2D[3002]  = tp_akns[3002].Thickness_2D
thick2D_t[3002] = thick2D[3002].sum(dim='time')
thick2D_timeseries_fix2 = np.array([])
for time_index in range(nt):        
        timeseries_value = np.sum(thick2D[3002][time_index,:,:].values)   
        thick2D_timeseries_fix2 = np.append(thick2D_timeseries_fix2, timeseries_value)

surfacemass_timeseries_fix2 = 500 * 440 *  mass_bal[3000][range(nt),i_density] * thick2D_timeseries_fix2 * 1e-9

surfacemass_timeseries_fix2.shape

plt.plot(surfacemass_timeseries_fix2[1:-1])
display(f'first value < 0 for nbrpartic[3000]: {thick2D_timeseries_fix2[0]} microns.  Plotting subsequent values.')
plt.ylabel('Integrated surface mass (tonnes)')
plt.xlabel('Time after oil release (hours) ')
plt.title('Spatially-integrated surface mass (tonnes) for run with biodegredation mass removed from MassOil')
plt.legend([f'3000-biofix (t(end) = {surfacemass_timeseries_fix2[nt-1]:3.1f})'])



StarcoderdataJupyterScriptsDedupFiltered