.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/example_droplet_distribution_compareJohansen2015.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_example_droplet_distribution_compareJohansen2015.py: Droplet distribution ================================== Plotting different droplet size distributions used in Opendrift .. GENERATED FROM PYTHON SOURCE LINES 8-15 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from datetime import datetime, timedelta from opendrift.models.openoil import OpenOil .. GENERATED FROM PYTHON SOURCE LINES 16-17 droplet size interval for plotting .. GENERATED FROM PYTHON SOURCE LINES 17-20 .. code-block:: Python dmin = 1.e-6 dmax = 3.e-3 .. GENERATED FROM PYTHON SOURCE LINES 21-22 simulation with Johansen et al. (2015) droplet spectrum .. GENERATED FROM PYTHON SOURCE LINES 22-47 .. code-block:: Python o = OpenOil(loglevel=20, weathering_model='noaa') o.set_config('environment:fallback:land_binary_mask', 0) o.set_config('environment:fallback:x_sea_water_velocity', -.2) o.set_config('environment:fallback:y_sea_water_velocity', 0) o.set_config('environment:fallback:x_wind', 6.) o.set_config('environment:fallback:y_wind', 0) o.set_config('environment:fallback:sea_water_temperature', 5.) o.set_config('environment:fallback:sea_surface_wave_stokes_drift_x_velocity', .3) o.set_config('environment:fallback:sea_surface_wave_stokes_drift_y_velocity', 0) o.set_config('wave_entrainment:droplet_size_distribution', 'Johansen et al. (2015)') o.set_config('processes:evaporation', False) o.set_config('processes:dispersion', False) o.seed_elements(lon=4, lat=60, time=datetime.utcnow(), number=10000, radius=100, z=0, oil_type='TROLL, STATOIL', oil_film_thickness=0.005) o.run(duration=timedelta(hours=2), time_step=3600) droplet_diameters = o.elements.diameter sd = 0.4 Sd = np.log(10.)*sd DV_50 = np.median(droplet_diameters) DN_50 = np.exp( np.log(DV_50) - 3*Sd**2 ) print('DV_50: %f' % DV_50) print('DN_50: %f' % DN_50) .. rst-class:: sphx-glr-script-out .. code-block:: none 13:15:12 INFO opendrift.models.basemodel:529: OpenDriftSimulation initialised (version 1.11.2 / v1.11.2-11-g21e4c60) 13:15:12 INFO opendrift.models.openoil.adios.dirjs:90: Querying ADIOS database for oil: TROLL, STATOIL 13:15:12 INFO opendrift.models.openoil.openoil:1721: Using density 891.412528 and viscosity 0.0032285361678330597 of oiltype TROLL, STATOIL 13:15:12 INFO opendrift.models.basemodel.environment:220: Adding a dynamical landmask with max. priority based on assumed maximum speed of 1.0 m/s. Adding a customised landmask may be faster... 13:15:17 INFO opendrift.models.basemodel.environment:247: Fallback values will be used for the following variables which have no readers: 13:15:17 INFO opendrift.models.basemodel.environment:250: x_sea_water_velocity: -0.200000 13:15:17 INFO opendrift.models.basemodel.environment:250: y_sea_water_velocity: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: x_wind: 6.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: y_wind: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_surface_height: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: upward_sea_water_velocity: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_surface_wave_significant_height: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_surface_wave_stokes_drift_x_velocity: 0.300000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_surface_wave_stokes_drift_y_velocity: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_surface_wave_period_at_variance_spectral_density_maximum: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_surface_wave_mean_period_from_variance_spectral_density_second_frequency_moment: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_ice_area_fraction: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_ice_x_velocity: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_ice_y_velocity: 0.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_water_temperature: 5.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_water_salinity: 34.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: sea_floor_depth_below_sea_level: 10000.000000 13:15:17 INFO opendrift.models.basemodel.environment:250: ocean_vertical_diffusivity: 0.020000 13:15:17 INFO opendrift.models.basemodel.environment:250: ocean_mixed_layer_thickness: 50.000000 13:15:17 INFO opendrift.models.basemodel:908: Using existing reader for land_binary_mask 13:15:17 INFO opendrift.models.basemodel:920: All points are in ocean 13:15:17 INFO opendrift.models.openoil.openoil:691: Oil-water surface tension is 0.032037 Nm 13:15:17 INFO opendrift.models.openoil.openoil:700: Using max water fractions [0.74, 0.75] for temperatures [5.0, 15.0] for oiltype TROLL, STATOIL 13:15:17 INFO opendrift.models.openoil.openoil:701: Corresponding max water fraction from GNOME is 0.5262214946654028 13:15:17 INFO opendrift.models.basemodel:2011: 2024-04-12 13:15:12.349239 - step 1 of 2 - 10000 active elements (0 deactivated) 13:15:18 INFO opendrift.models.basemodel:2011: 2024-04-12 14:15:12.349239 - step 2 of 2 - 10000 active elements (0 deactivated) /root/project/examples/example_droplet_distribution_compareJohansen2015.py:42: RuntimeWarning: divide by zero encountered in log DN_50 = np.exp( np.log(DV_50) - 3*Sd**2 ) DV_50: 0.000000 DN_50: 0.000000 .. GENERATED FROM PYTHON SOURCE LINES 48-49 Plotting .. GENERATED FROM PYTHON SOURCE LINES 49-62 .. code-block:: Python plt.figure(figsize=[14,14]) plt.subplot(3, 2, 1) nVpdf, binsV, patches = plt.hist(droplet_diameters, 100, range=(dmin,dmax), align='mid') plt.xlabel('Droplet diameter d [m]', fontsize=8) plt.ylabel('V(d)', fontsize=8) plt.title('volume spectrum\nJohansen et al. (2015) distribution in OpenOil', fontsize=10) plt.subplot(3,2,2) nVcum, binsV, patches = plt.hist(droplet_diameters, 100, range=(dmin,dmax), align='mid', cumulative=True) plt.xlabel('Droplet diameter d [m]', fontsize=8) plt.ylabel('V(d)', fontsize=8) plt.title('cumulative volume spectrum', fontsize=10) .. image-sg:: /gallery/images/sphx_glr_example_droplet_distribution_compareJohansen2015_001.png :alt: volume spectrum Johansen et al. (2015) distribution in OpenOil, cumulative volume spectrum :srcset: /gallery/images/sphx_glr_example_droplet_distribution_compareJohansen2015_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'cumulative volume spectrum') .. GENERATED FROM PYTHON SOURCE LINES 63-64 calculate number spectrum from volume spectrum .. GENERATED FROM PYTHON SOURCE LINES 64-69 .. code-block:: Python d = 0.5* (binsV[1:] + binsV[:-1]) V = 4./3. * np.pi * (d/2.)**3 Npdf = nVpdf / V Ncum = np.cumsum(Npdf) .. GENERATED FROM PYTHON SOURCE LINES 70-71 calculate theoretical cumulative Number distribution .. GENERATED FROM PYTHON SOURCE LINES 71-119 .. code-block:: Python spectrum = (np.exp(-(np.log(d) - np.log(DN_50))**2 / (2 * Sd**2))) / (d * Sd * np.sqrt(2 * np.pi)) # from OpenOil median diameter DN_50_johansen = 0.000483 spectrumJ = (np.exp(-(np.log(d) - np.log(DN_50_johansen))**2 / (2 * Sd**2))) / (d * Sd * np.sqrt(2 * np.pi)) # from parameters in Johansen et al. 2015 spectrum_cum = np.cumsum(spectrum) spectrumJ_cum = np.cumsum(spectrumJ) # calculate theoretical number distribution pdf #spectrum_pdf = spectrum / np.sum(spectrum) #spectrumJ_pdf = spectrumJ / np.sum(spectrumJ) plt.subplot(3,2,3) plt.plot(d, Npdf/np.sum(Npdf), lw=2) #plt.plot(d, spectrum/np.sum(spectrum), label='curve fit from OpenOil', lw=2) plt.plot(d, spectrumJ/np.sum(spectrumJ), label='curve fit from Johansen et al. 2015', lw=2) plt.xlabel('Droplet diameter d [m]', fontsize=8) plt.ylabel('N(d)', fontsize=8) plt.title('number spectrum', fontsize=10) plt.subplot(3,2,4) plt.plot(d, Ncum/np.max(Ncum), label='OpenOil result', lw=2) #plt.plot(d, spectrum_cum/np.max(spectrum_cum), label='OpenOil par.', lw=2) plt.plot(d, spectrumJ_cum/np.max(spectrumJ_cum), label='Johansen et al. 2015', lw=2) plt.xlabel('Droplet diameter d [m]', fontsize=8) plt.ylabel('N(d)', fontsize=8) plt.title('cumulative number spectrum', fontsize=10) plt.legend(loc='lower right') plt.subplot(3,2,5) plt.loglog(d, Npdf/np.sum(Npdf), lw=2) #plt.loglog(d, spectrum/np.sum(spectrum), label='curve fit from OpenOil', lw=2) plt.loglog(d, spectrumJ/np.sum(spectrumJ), label='curve fit from Johansen et al. 2015', lw=2) plt.xlabel('Droplet diameter d [m]', fontsize=8) plt.ylabel('N(d)', fontsize=8) plt.title('number spectrum', fontsize=10) plt.subplot(3,2,6) plt.semilogx(d, Ncum/np.max(Ncum), label='OpenOil result', lw=2) #plt.semilogx(d, spectrum_cum/np.max(spectrum_cum), label='OpenOil par.', lw=2) plt.semilogx(d, spectrumJ_cum/np.max(spectrumJ_cum), label='Johansen et al. 2015', lw=2) plt.xlabel('Droplet diameter d [m]', fontsize=8) plt.ylabel('N(d)', fontsize=8) plt.title('cumulative number spectrum', fontsize=10) #plt.legend(loc='lower right') plt.show() .. image-sg:: /gallery/images/sphx_glr_example_droplet_distribution_compareJohansen2015_002.png :alt: number spectrum, cumulative number spectrum, number spectrum, cumulative number spectrum :srcset: /gallery/images/sphx_glr_example_droplet_distribution_compareJohansen2015_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /root/project/examples/example_droplet_distribution_compareJohansen2015.py:71: RuntimeWarning: divide by zero encountered in log spectrum = (np.exp(-(np.log(d) - np.log(DN_50))**2 / (2 * Sd**2))) / (d * Sd * np.sqrt(2 * np.pi)) # from OpenOil median diameter .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 8.314 seconds) .. _sphx_glr_download_gallery_example_droplet_distribution_compareJohansen2015.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_droplet_distribution_compareJohansen2015.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_droplet_distribution_compareJohansen2015.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_