_GenerateBCRecordsOldUgly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write discharge records\n",
    "### ⚠️ Writes only records part (for two sediments) of `.bct` or `.bcc` files! ⚠️ \n",
    "\n",
    "#### Required input \n",
    "* timestep\n",
    "* start and end of discharge time\n",
    "* Number of sigma layers\n",
    "* Number of sigma layer over which boundary condition is imposed\n",
    "* Sediment weight per volume [kg/m3]\n",
    "\n",
    "\n",
    "## TODO\n",
    "* Split into a function\n",
    "    * Multiple sediment + multiple times!\n",
    "* Use pandas\n",
    "* Don't hardcode end time as string\n",
    "* Read vertical grid from .mdf and .dep file\n",
    "* Read horizontal cell size from grid file\n",
    "* Plot time vs discharge of boundary condition\n",
    "* Make more modular to allow for more than two sediments\n",
    "    * Make a function for writing records to new file\n",
    "* Add some plotting of the column to visualize height of layers and among which layers flow is imposed\n",
    "* Write the whole .bcc/.bct files instead of just the records\n",
    "* Read values from mdf and dep file\n",
    "\n",
    "## Might be good\n",
    "1. Use [Pint](https://pint.readthedocs.io/en/0.11/) to keep track of units\n",
    "```python\n",
    "import pint\n",
    "ureg = pint.UnitRegistry()\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from JulesD3D.utils import formatSci\n",
    "from IPython.display import Markdown as md"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nr_sigma_layers = 80 # total sigma layers\n",
    "nr_discharge_layers = 53\n",
    "nr_zero_discharge_layers = nr_sigma_layers - nr_discharge_layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sand_mass_conc = 19.875 # [kg/m3] #  # 25% 39.74\n",
    "silt_mass_conc = 59.625 # [kg/m3] #  # 75% 39.74"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_width = 200   # ureg.meter [m] \n",
    "cell_length = 200  # ureg.meter [m]\n",
    "init_depth = 300   # ureg.meter [m]\n",
    "density_sed = 2650 # [kg/m3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "discharge_total = 4500 # [m3/s]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.2300000e+003'"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "timestep = 0.3 # [min]\n",
    "start_discharge_time = 9 # [min]\n",
    "discharge_duration = 39 # [min]\n",
    "\n",
    "end_time = 1230\n",
    "formatSci(end_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "nr_of_timesteps = int(end_time / timestep)\n",
    "nr_of_timesteps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "39"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "discharge_duration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "end_discharge_time = start_discharge_time + discharge_duration # [min]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Discharge is linearly interpolated between timesteps\n",
    "Therefore, we need to add one timestep before the start and one timestep after in which discharges are zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pre_start_discharge_time = start_discharge_time - timestep  # [min]\n",
    "post_end_discharge_time = end_discharge_time + timestep  # [min]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| ⚠️ WARNING: Always check wether time is a multiple of timestep! ⚠️ |\n",
    "| --- |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "\n",
       "# Sanity check\n",
       "Weird check for floating point errors\n",
       "## Is discharge duration a multiple of timestep? True!\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from decimal import Decimal\n",
    "md(f'''\n",
    "# Sanity check\n",
    "Weird check for floating point errors\n",
    "## Is discharge duration a multiple of timestep? {Decimal(discharge_duration) % Decimal(timestep) < 10e-9}!\n",
    "''') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# should be read from mdf file\n",
    "sigma_layers_percentages = np.array([7.180, 6.000, 5.500, 5.000, 4.500, 4.000, 3.800, 3.600, 3.400, 3.200, 3.000, 2.900, 2.800, 2.700, 2.600, 2.500, 2.400, 2.300, 2.200, 2.100, 2.000, 1.900, 1.800, 1.700, 1.600, 1.500, 1.400, 1.300, 1.200, 1.100, 1.000, 0.950, 0.900, 0.850, 0.800, 0.750, 0.700, 0.650, 0.600, 0.550, 0.500, 0.450, 0.400, 0.350, 0.300, 0.270, 0.230, 0.200, 0.150, 0.130, 0.120, 0.110, 0.100, 0.100, 0.090, 0.090, 0.090, 0.090, 0.090, 0.080, 0.080, 0.080, 0.080, 0.070, 0.070, 0.070, 0.060, 0.060, 0.060, 0.060, 0.060, 0.050, 0.050, 0.050, 0.050, 0.040, 0.040, 0.040, 0.030, 0.03])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['0.0000000e+000',\n",
       " '8.7000000e+000',\n",
       " '9.0000000e+000',\n",
       " '4.8000000e+001',\n",
       " '4.8300000e+001',\n",
       " '1.2300000e+003']"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# TODO: hardcoded simultion start + edn time string\n",
    "times = [f\"{formatSci(0)}\", \n",
    "         f\"{formatSci(pre_start_discharge_time)}\", \n",
    "         f\"{formatSci(start_discharge_time)}\",\\\n",
    "         f\"{formatSci(end_discharge_time)}\", \n",
    "         f\"{formatSci(post_end_discharge_time)}\", \n",
    "         formatSci(end_time)] # \"1.2300000e+003\"]\n",
    "times"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_zeros_line = '  '.join('0.0000000e+000' for i in range(nr_sigma_layers))\n",
    "line_all_zero = f\"{all_zeros_line}  {all_zeros_line}\" # A + B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repeat '0.0000000e+000' nr_zero_discharge_layers times\n",
    "# np.zeroes(zero_records_discharge_layer)\n",
    "zero_records_discharge_layer = '  '.join('0.0000000e+000' for i in range(nr_zero_discharge_layers))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write Concentration file (.bcc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "formatted_sand_conc = formatSci(sand_mass_conc)\n",
    "formatted_silt_conc = formatSci(silt_mass_conc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "### The first 27 'cells' are all zeroes, then the next 53 layers contain the discharge values"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "md(f\"### The first {nr_zero_discharge_layers} 'cells' are all zeroes, then the next {nr_discharge_layers} layers contain the discharge values\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: should loop over arrray of sediments and then write each line but for only two like now it works fine\n",
    "sand_conc_records_discharge_layer = '  '.join(formatted_sand_conc for i in range(nr_discharge_layers))\n",
    "silt_conc_records_discharge_layer = '  '.join(formatted_silt_conc for i in range(nr_discharge_layers))\n",
    "\n",
    "discharge_sand_vol_line = f\"{zero_records_discharge_layer}  {sand_conc_records_discharge_layer}\"\n",
    "discharge_silt_vol_line = f\"{zero_records_discharge_layer}  {silt_conc_records_discharge_layer}\"\n",
    "\n",
    "# Each line has two discharge 'points'\n",
    "line_with_sand_conc_records = f'{discharge_sand_vol_line}  {discharge_sand_vol_line}' # A + B\n",
    "line_with_silt_conc_records = f'{discharge_silt_vol_line}  {discharge_silt_vol_line}' # A + B\n",
    "bcc_filename = './generated/7525_15.bccrecords'\n",
    "\n",
    "# time + A + B\n",
    "with open(bcc_filename, 'w') as f:\n",
    "    # why not use a f multiline string here f'''{ line below} '''\n",
    "    linesToWrite = [\n",
    "        \"sand\\n\",\n",
    "        f\" {times[0]}  {line_all_zero}\" + '\\n',\n",
    "        f\" {times[1]}  {line_all_zero}\" + '\\n',\n",
    "        f\" {times[2]}  {line_with_sand_conc_records}\" + '\\n',\n",
    "        f\" {times[3]}  {line_with_sand_conc_records}\" + '\\n',\n",
    "        f\" {times[4]}  {line_all_zero}\" + '\\n',\n",
    "        f\" {times[5]}  {line_all_zero}\" + '\\n',\n",
    "        \"\\n\",\n",
    "        \"silt\\n\"\n",
    "        f\" {times[0]}  {line_all_zero}\" + '\\n',\n",
    "        f\" {times[1]}  {line_all_zero}\" + '\\n',\n",
    "        f\" {times[2]}  {line_with_silt_conc_records}\" + '\\n',\n",
    "        f\" {times[3]}  {line_with_silt_conc_records}\" + '\\n',\n",
    "        f\" {times[4]}  {line_all_zero}\" + '\\n',\n",
    "        f\" {times[5]}  {line_all_zero}\" + '\\n'\n",
    "    ]\n",
    "    f.writelines(line for line in linesToWrite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write Discharge file (.bct)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "discharge_layer_sigma_percentages = sigma_layers_percentages[nr_zero_discharge_layers:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Discharge layers: 16.42%"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "discharge_layers_sigma_percentage_sum = discharge_layer_sigma_percentages.sum()\n",
    "md(f\"Discharge layers: {discharge_layers_sigma_percentage_sum}%\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Now divide 4500 $m^3/s$ discharge among the bottom 16.42% of sigma layers"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "md(f\"Now divide {discharge_total} $m^3/s$ discharge among the bottom {discharge_layers_sigma_percentage_sum}% of sigma layers\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_discharge_layers = discharge_layer_sigma_percentages/discharge_layers_sigma_percentage_sum * discharge_total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Check if its right:\n",
       "    Sum of discharges over 53 discharge layers is 4499.999999999998 $m^3/s$\n",
       "    and entered discharge is 4500 $m^3/s$, so it's **not right**\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "correct = discharge_total == new_discharge_layers.sum()\n",
    "\n",
    "md(f'''Check if its right:\n",
    "    Sum of discharges over {nr_discharge_layers} discharge layers is {new_discharge_layers.sum()} $m^3/s$\n",
    "    and entered discharge is {discharge_total} $m^3/s$, so it's {\"**correct!**\" if correct else \"**not right man**\"}\n",
    "    ''')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "discharge_records_discharge_layer = '  '.join(formatSci(discharge) for discharge in new_discharge_layers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "discharge_line_formatted = f\"{zero_records_discharge_layer}  {discharge_records_discharge_layer}\"\n",
    "# Each line has two discharge 'points'\n",
    "line_with_records = f'{discharge_line_formatted}  {discharge_line_formatted}' # A + B\n",
    "bct_filename = f'./generated/discharge{discharge_total}.bctrecords'\n",
    "\n",
    "# Format: time + A + B\n",
    "with open(bct_filename, 'w') as f:\n",
    "    linesToWrite = [\n",
    "        f\" {times[0]}  {line_all_zero}\\n\",\n",
    "        f\" {times[1]}  {line_all_zero}\\n\",\n",
    "        f\" {times[2]}  {line_with_records}\\n\",\n",
    "        f\" {times[3]}  {line_with_records}\\n\",\n",
    "        f\" {times[4]}  {line_all_zero}\\n\",\n",
    "        f\" {times[5]}  {line_all_zero}\\n\"\n",
    "    ]\n",
    "    f.writelines(line for line in linesToWrite)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# should have concat'ed the values like this first then formatting instead of joining formatted stuff\n",
    "discharge_line_values = np.hstack([np.zeros(nr_zero_discharge_layers), new_discharge_layers]) # m3/s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate a bunch-o-stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_heights = np.array(sigma_layers_percentages)/100 * init_depth # [m]\n",
    "cell_volumes = cell_heights * cell_width * cell_length # [m3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate concentrations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate total sand volume in column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "sand_zero_dis_layer = np.zeros(nr_zero_discharge_layers) # [kg/m3]\n",
    "sand_dis_layer = np.ones(nr_discharge_layers) * sand_mass_conc # [kg/m3]\n",
    "\n",
    "all_sand_layers_weight_vol = np.hstack([sand_zero_dis_layer, sand_dis_layer]) # [kg/m3]\n",
    "vol_sand_per_layer = all_sand_layers_weight_vol * cell_volumes / density_sed # [m3]\n",
    "total_sand_vol_in_column = vol_sand_per_layer.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate total silt volume in column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "silt_zero_dis_layer = np.zeros(nr_zero_discharge_layers)\n",
    "silt_dis_layer = np.ones(nr_discharge_layers) * silt_mass_conc\n",
    "\n",
    "all_silt_layers_weight_vol = np.hstack([silt_zero_dis_layer, silt_dis_layer]) # [kg/m3]\n",
    "vol_silt_per_layer = all_silt_layers_weight_vol * cell_volumes / density_sed # [m3]\n",
    "total_silt_vol_in_column = vol_silt_per_layer.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_vol_discharge_cells = cell_volumes[nr_zero_discharge_layers:].sum() # [m3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "sand_conc = total_sand_vol_in_column/total_vol_discharge_cells\n",
    "silt_conc = total_silt_vol_in_column/total_vol_discharge_cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_sediment_conc = (silt_conc+sand_conc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate total sediment weight and volume coming in\n",
    "$[m^3/s] \\times [kg/m^3] \\times [s] = [kg]$\n",
    "\n",
    "discharge $\\times$ weight volume $\\times$ duration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "sand_weight_flowing_in = discharge_line_values * all_sand_layers_weight_vol * discharge_duration * 60 # LIST [m3/s] * [kg/m3] * [s] = [kg] \n",
    "total_sand_weight_flowing_in = sand_weight_flowing_in.sum() # [kg]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "silt_weight_flowing_in = discharge_line_values * all_silt_layers_weight_vol * discharge_duration * 60 # LIST [m3/s] * [kg/m3] * [s] = [kg] \n",
    "total_silt_weight_flowing_in = silt_weight_flowing_in.sum() # [kg]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_sediment_weight_flowing_in = total_sand_weight_flowing_in + total_silt_weight_flowing_in\n",
    "total_sediment_volume_flowing_in = total_sediment_weight_flowing_in/density_sed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "\n",
       "* **Total sediment volume flowing in = 3.159 $\\times10^5 m^3$**\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "md(f'''\n",
    "* **Total sediment volume flowing in = {total_sediment_volume_flowing_in/10**5} $\\\\times10^5 m^3$**\n",
    "''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_discharge_volume = discharge_duration * 60 * discharge_total # [s * m3/s = m3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "\n",
       "# Summary\n",
       "\n",
       "### Time\n",
       "* Total simulated time = 20.5 hours\n",
       "* Number of timesteps = 4100\n",
       "\n",
       "\n",
       "\n",
       "### Flow\n",
       "* Total flow discharge volume = 10.53 $ \\times 10^6 m^3$\n",
       "* Discharge height = 49.26 m\n",
       "* Total discharge = 4500 $m^3/s$\n",
       "\n",
       "* Discharge start time = 9\n",
       "* Discharge end time = 48\n",
       "* Discharge duration = 39 min\n",
       "\n",
       "### Sediments\n",
       "* Silt weight per volume = 59.625 $kg/m^3$\n",
       "* Sand weight per volume = 19.875 $kg/m^3$\n",
       "\n",
       "* Total silt volume 0.79 $\\times10^5 m^3$\n",
       "* Total sand volume 2.37 $\\times10^5 m^3$\n",
       "* **Total sediment volume = 3.159 $\\times10^5 m^3$**\n",
       "\n",
       "* Silt Concentration = 2.25%\n",
       "* Sand Concentration =  0.75%\n",
       "* Total sediment concentration = 3.0%\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# turn this into a markdown table\n",
    "md(f'''\n",
    "# Summary\n",
    "\n",
    "### Time\n",
    "* Total simulated time = {int(end_time)/60} hours\n",
    "* Number of timesteps = {nr_of_timesteps}\n",
    "\n",
    "\n",
    "\n",
    "### Flow\n",
    "* Total flow discharge volume = {total_discharge_volume/10**6} $ \\\\times 10^6 m^3$\n",
    "* Discharge height = {cell_heights[nr_zero_discharge_layers:].sum()} m\n",
    "* Total discharge = {discharge_total} $m^3/s$\n",
    "\n",
    "* Discharge start time = {start_discharge_time}\n",
    "* Discharge end time = {end_discharge_time}\n",
    "* Discharge duration = {end_discharge_time - start_discharge_time} min\n",
    "\n",
    "### Sediments\n",
    "* Silt weight per volume = {silt_mass_conc} $kg/m^3$\n",
    "* Sand weight per volume = {sand_mass_conc} $kg/m^3$\n",
    "\n",
    "* Total silt volume {round(total_sand_weight_flowing_in/density_sed/10**5, 2)} $\\\\times10^5 m^3$\n",
    "* Total sand volume {round(total_silt_weight_flowing_in/density_sed/10**5, 2)} $\\\\times10^5 m^3$\n",
    "* **Total sediment volume = {round(total_sediment_volume_flowing_in/10**5, 3)} $\\\\times10^5 m^3$**\n",
    "\n",
    "* Silt Concentration = {round(silt_conc * 100, 2)}%\n",
    "* Sand Concentration =  {round(sand_conc * 100, 2)}%\n",
    "* Total sediment concentration = {round(total_sediment_conc * 100, 2)}%\n",
    "\n",
    "''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<__main__.MonsterBCCClass at 0x1145eb210>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# shouldnt it get the sigma stuff from an mdf class or something?\n",
    "class MonsterBCCClass():\n",
    "    def __init__(self, nr_sigma_layers, nr_discharge_layers):\n",
    "        self.nr_sigma_layers = nr_sigma_layers\n",
    "        self.nr_discharge_layers = nr_discharge_layers\n",
    "        self.nr_zero_discharge_layers = nr_sigma_layers - nr_discharge_layers\n",
    "        # self.sigma_percentage = sigma_percentage\n",
    "        \n",
    "        print(self.nr_zero_discharge_layers)\n",
    "MonsterBCCClass(80, 57)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # mass conc is just repeated nr_discharge_layers times\n",
    "# # whereas discharge values are cell specific \n",
    "# def makeDischargeRecords(mass_conc, nr_sigma_layers=80, nr_discharge_layers = 53):\n",
    "#     # remove to class\n",
    "#     nr_zero_discharge_layers = nr_sigma_layers - nr_discharge_layers    \n",
    "#     formatted_mass_conc = formatSci(mass_conc)\n",
    "#     #####\n",
    "#     print(f\"The first {nr_zero_discharge_layers} 'cells' are all zeroes, then the next {nr_discharge_layers} layers contain the discharge values: {mass_conc}\")\n",
    "    \n",
    "#     zero_records_discharge_layer = '  '.join('0.0000000e+000' for i in range(nr_zero_discharge_layers))\n",
    "    \n",
    "#     # what if discharge layers are an argument to this function?\n",
    "#     records_discharge_layers = '  '.join(formatted_mass_conc for i in range(nr_discharge_layers))\n",
    "    \n",
    "#     # first the cells where discharge is zero, then the cells that discharge \n",
    "#     discharge_line = f\"{zero_records_discharge_layer}  {records_discharge_layers}\"\n",
    "\n",
    "#     # Each line has two discharge 'points' so repeat it twice\n",
    "#     complete_line_conc_records = f'{discharge_line}  {discharge_line}' # A + B\n",
    "\n",
    "#     return complete_line_conc_records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def allZeroesRecords(nr_sigma_layers=80):\n",
    "    all_zeros_line = '  '.join('0.0000000e+000' for i in range(nr_sigma_layers))\n",
    "    line_all_zero = f\"{all_zeros_line}  {all_zeros_line}\" # left + right\n",
    "    \n",
    "    # Repeat '0.0000000e+000' nr_zero_discharge_layers times\n",
    "    zero_records_discharge_layer = '  '.join('0.0000000e+000' for i in range(nr_zero_discharge_layers))\n",
    "    \n",
    "    return zero_records_discharge_layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeBcTimes(timestep=0.3, start_discharge_time=15, discharge_duration=45):\n",
    "    # Discharge is linearly interpolated between timesteps\n",
    "    # Therefore, we need to add one timestep before the start and one timestep after in which discharges are zero\n",
    "    end_discharge_time = start_discharge_time + discharge_duration # [min]\n",
    "    pre_start_discharge_time = start_discharge_time - timestep  # [min]\n",
    "    post_end_discharge_time = end_discharge_time + timestep  # [min]    \n",
    "\n",
    "    times = [formatSci(0), \n",
    "             formatSci(pre_start_discharge_time), \n",
    "             formatSci(start_discharge_time),\n",
    "             formatSci(end_discharge_time), \n",
    "             formatSci(post_end_discharge_time), \n",
    "             formatSci(end_time)\n",
    "            ]\n",
    "    \n",
    "    return times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcDischargePerCell(discharge_total=4500, sigma_layers_percentages=[], nr_zero_discharge_layers=0):\n",
    "    discharge_layer_sigma_percentages = sigma_layers_percentages[nr_zero_discharge_layers:] # select only the discharge parts\n",
    "    discharge_layers_sigma_percentage_sum = discharge_layer_sigma_percentages.sum() # sum it\n",
    "\n",
    "    print(f\"Now divide {discharge_total} $m^3/s$ discharge among the bottom {discharge_layers_sigma_percentage_sum}% height of sigma layers\")\n",
    "    new_discharge_layers = discharge_layer_sigma_percentages/discharge_layers_sigma_percentage_sum * discharge_total\n",
    "    \n",
    "    discharge_records = '  '.join(formatSci(discharge) for discharge in new_discharge_layers)\n",
    "\n",
    "    return discharge_records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeDischargeRecords(discharge_records=[], nr_sigma_layers=80, nr_discharge_layers = 53):\n",
    "    # remove to class\n",
    "    nr_zero_discharge_layers = nr_sigma_layers - nr_discharge_layers    \n",
    "\n",
    "    print(f\"The first {nr_zero_discharge_layers} 'cells' are all zeroes, then the next {nr_discharge_layers} layers contain the discharge values\")\n",
    "    \n",
    "    zero_records_discharge_layer = '  '.join('0.0000000e+000' for i in range(nr_zero_discharge_layers))\n",
    "    \n",
    "    # first the cells where discharge is zero, then the cells that discharge \n",
    "    discharge_line = f\"{zero_records_discharge_layer}  {discharge_records}\"\n",
    "\n",
    "    # Each line has two discharge 'points' so repeat it twice\n",
    "    complete_line_conc_records = f'{discharge_line}  {discharge_line}' # A + B\n",
    "\n",
    "    return complete_line_conc_records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now divide 4500 $m^3/s$ discharge among the bottom 16.42% height of sigma layers\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'3.5627284e+002  3.2886724e+002  3.0146163e+002  2.7405603e+002  2.6035323e+002  2.4665043e+002  2.3294762e+002  2.1924482e+002  2.0554202e+002  1.9183922e+002  1.7813642e+002  1.6443362e+002  1.5073082e+002  1.3702801e+002  1.2332521e+002  1.0962241e+002  9.5919610e+001  8.2216809e+001  7.3995128e+001  6.3032887e+001  5.4811206e+001  4.1108404e+001  3.5627284e+001  3.2886724e+001  3.0146163e+001  2.7405603e+001  2.7405603e+001  2.4665043e+001  2.4665043e+001  2.4665043e+001  2.4665043e+001  2.4665043e+001  2.1924482e+001  2.1924482e+001  2.1924482e+001  2.1924482e+001  1.9183922e+001  1.9183922e+001  1.9183922e+001  1.6443362e+001  1.6443362e+001  1.6443362e+001  1.6443362e+001  1.6443362e+001  1.3702801e+001  1.3702801e+001  1.3702801e+001  1.3702801e+001  1.0962241e+001  1.0962241e+001  1.0962241e+001  8.2216809e+000  8.2216809e+000'"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calcDischargePerCell(sigma_layers_percentages, nr_zero_discharge_layers=nr_zero_discharge_layers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The first 27 'cells' are all zeroes, then the next 53 layers contain the discharge values\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  0.0000000e+000  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001  1.9875000e+001'"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formatted_mass_conc = formatSci(sand_mass_conc)\n",
    "conc_records_discharge_layers = '  '.join(formatted_mass_conc for i in range(nr_discharge_layers))\n",
    "makeDischargeRecords(conc_records_discharge_layers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc_filename = './generated/7525_15.bccrecords'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "non-default argument follows default argument (<ipython-input-19-97e09f30812f>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-19-97e09f30812f>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    def generateVerticalBC(timestep = 0.3, nr_sigma_layers=80, start_discharge_time, end_discharge_time, sand_mass_conc=19.75, silt_mass_conc=19.75):\u001b[0m\n\u001b[0m                          ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m non-default argument follows default argument\n"
     ]
    }
   ],
   "source": [
    "def generateVerticalBC(filename='untitled.bccrecords', timestep = 0.3, nr_sigma_layers=80, start_discharge_time=0, discharge_duration=0, sand_mass_conc=19.75, silt_mass_conc=19.75):\n",
    "    if discharge_duration == 0:\n",
    "        raise Exception(\"discharge_duration cannot be 0\")\n",
    "    \n",
    "    if not Decimal(discharge_duration) % Decimal(timestep) < 10e-9:\n",
    "        raise Exception(\"Discharge duration is not a multiple of timestep\")\n",
    "\n",
    "    times = makeBcTimes(timestep=timestep, start_discharge_time=start_discharge_time, end_discharge_time=end_discharge_time)\n",
    "\n",
    "    line_with_sand_conc_records = makeDischargeRecords(sand_mass_conc, nr_sigma_layers, nr_discharge_layers)\n",
    "    line_with_silt_conc_records = makeDischargeRecords(silt_mass_conc, nr_sigma_layers, nr_discharge_layers)\n",
    "    \n",
    "    # --- Write Concentration file (.bcc) ---\n",
    "    # time + A + B\n",
    "    with open(bcc_filename, 'w') as f:\n",
    "        # why not use a f multiline string here f'''{ line below} '''\n",
    "        # what if more than 6 records? this approach is inflexible\n",
    "        # for sediment in sediments:\n",
    "            # for time in times:\n",
    "            \n",
    "            # one discharge always follows this sequence of 6 time values with discharge records at the 3rd and 4th place\n",
    "            \n",
    "        linesToWrite = [\n",
    "            \"sand\\n\",\n",
    "            f\" {times[0]}  {line_all_zero}\" + '\\n',\n",
    "            f\" {times[1]}  {line_all_zero}\" + '\\n',\n",
    "            f\" {times[2]}  {line_with_sand_conc_records}\" + '\\n',\n",
    "            f\" {times[3]}  {line_with_sand_conc_records}\" + '\\n',\n",
    "            f\" {times[4]}  {line_all_zero}\" + '\\n',\n",
    "            f\" {times[5]}  {line_all_zero}\" + '\\n',\n",
    "            \"\\n\",\n",
    "            \"silt\\n\"\n",
    "            f\" {times[0]}  {line_all_zero}\" + '\\n',\n",
    "            f\" {times[1]}  {line_all_zero}\" + '\\n',\n",
    "            f\" {times[2]}  {line_with_silt_conc_records}\" + '\\n',\n",
    "            f\" {times[3]}  {line_with_silt_conc_records}\" + '\\n',\n",
    "            f\" {times[4]}  {line_all_zero}\" + '\\n',\n",
    "            f\" {times[5]}  {line_all_zero}\" + '\\n'\n",
    "        ]\n",
    "        f.writelines(line for line in linesToWrite)\n",
    "        \n",
    "    bct_filename = f'./generated/discharge{discharge_total}.bctrecords'\n",
    "\n",
    "    # Format: time + A + B\n",
    "    with open(bct_filename, 'w') as bct_file:\n",
    "        linesToWrite = [\n",
    "            f\" {times[0]}  {line_all_zero}\\n\",\n",
    "            f\" {times[1]}  {line_all_zero}\\n\",\n",
    "            f\" {times[2]}  {line_with_records}\\n\",\n",
    "            f\" {times[3]}  {line_with_records}\\n\",\n",
    "            f\" {times[4]}  {line_all_zero}\\n\",\n",
    "            f\" {times[5]}  {line_all_zero}\\n\"\n",
    "        ]\n",
    "        bct_file.writelines(line for line in linesToWrite) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}