Cortix 2019 07Aug2019

Droplet Swirl Example

  • This is part of the Cortix-on-Jupyter-Notebook examples.
  • You must be in a Jupyter Notebook server to run this notebook.
  • Select each of the cells below and run them sequentially (use the run button, >| on the tool bar or use the Cell option on the menu bar).
  • Alternatively, on the menu bar run all cells: Cell -> Run All.

$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\vvar}{\boldsymbol{v}} \newcommand{\fvar}{\boldsymbol{f}} \newcommand{\Power}{\mathcal{P}} \newcommand{\bm}[1]{{\boldsymbol{#1}}} $

Introduction

This Cortix use-case simulates the motion of a swarm of droplets in a vortex stream. It consists of two modules, namely, a Droplet module used to model the droplet dynamics, and a Vortex module used to model the effects of the surrounding air on the falling droplets. The Droplet module is instantiated as many times as there are droplets in the simulation while a single Vortex module is connected to all Droplet instances. The communication between modules entails a two-way data exchange between the Vortex module and the Droplet modules, where Droplet sends
its position to Vortex and Vortex returns the air velocity to Droplet at the given position.

Droplet motion model

The equation of motion of a spherical droplet can be written as: \begin{equation*} m_\text{d}\,d_t\bm{v} = \bm{f}_\text{d} + \bm{f}_\text{b} , \end{equation*} where \begin{equation*} \bm{f}_\text{d} = c_\text{d} A \, \rho_\text{f}\, \frac{||\bm{v} - \bm{v}_\text{f}||}{2}\,(\bm{v} - \bm{v}_\text{f}) , \end{equation*} is the form drag force on the droplet, \begin{equation*} \bm{f}_\text{b} = (m_\text{d} - m_\text{f})\,g \hat{z} , \end{equation*} is the buoyancy force on the droplet, \begin{equation*} c_\text{d}(Re) = \begin{cases} \frac{24}{Re} & Re < 0.1\\ \Bigl(\sqrt{\frac{24}{Re}} + 0.5407\Bigr)^2 & 0.1 \leq Re < 6000 \\ 0.44 & Re \geq 6000 \end{cases} \end{equation*} is the drag coefficient as a function of Reynold's number, $Re=\frac{\rho_\text{f}\,||\bm{v}||\,d}{\mu_\text{f}}$. The mass of the droplet and its displaced fluid mass are denoted $m_d$ and $m_f$, respectively. Droplet diameter, $d$, dynamic viscosity, $\mu_\text{f}$, and mass density, $\rho_\text{f}$, of the surrounding air are provided.

Vortex motion model

Here we simply use an imposed vortex circulation in analytical form given by its tangential component of velocity \begin{equation*} v_\theta(r,z,t) = \Bigl(1 - e^{\frac{-r^2}{8\,r_c^2}}\Bigr) \frac{\Gamma}{2\pi\, \max(r,r_c)} f(z) \, \bigl|\cos(\mu\,t)\bigr|, \end{equation*} and its vertical component \begin{equation*} v_z(z,t) = v_h \, f(z) \, \bigl|\cos(\mu\,t)\bigr|, \end{equation*} where \begin{equation*} f(z) = e^{\frac{-(h - z)}{\ell}} \end{equation*} is a vertical relaxation factor, $r_c$ is the vortex core radius, $\Gamma = \frac{2\pi R}{v_\theta |_{r = R}}$ is the vortex circulation, $R$ is the vortex outer radius, $h$ is the height of the vortex, and $\ell$ is the relaxation length of $v_z$.

Write the run context

In [1]:
# Import various packages; must have the Cortix repository installed

import scipy.constants as const
import matplotlib.pyplot as plt

# Leave this block here for Azure
try:
    import cortix
except ImportError:
    print('Installing the "cortix" package...')
    print('')
    !pip install cortix
    import cortix

from cortix.src.cortix_main import Cortix
from cortix.src.network import Network

# Import the example modules
from cortix.examples.droplet_swirl.vortex import Vortex
from cortix.examples.droplet_swirl.droplet import Droplet
In [2]:
# Create a Cortix object with Python multiprocessing

swirl = Cortix(use_mpi=False,splash=True)

swirl.network = Network()
swirl_net = swirl.network

# Set parameters in SI units

n_droplets = 15
end_time = 3*const.minute
time_step = 0.2
[65266] 2019-08-28 09:54:34,471 - cortix - INFO - Created Cortix object 
_____________________________________________________________________________
                             L A U N C H I N G                               
_____________________________________________________________________________
      ...                                        s       .     (TAAG Fraktur)
   xH88"`~ .x8X                                 :8      @88>
 :8888   .f"8888Hf        u.      .u    .      .88      %8P      uL   ..
:8888>  X8L  ^""`   ...ue888b   .d88B :@8c    :888ooo    .     [email protected]  @88R
X8888  X888h        888R Y888r ="8888f8888r -*8888888  [email protected]  ""Y888k/"*P
88888  !88888.      888R I888>   4888>"88"    8888    888E`    Y888L
88888   %88888      888R I888>   4888> "      8888      888E      8888
88888 `> `8888>     888R I888>   4888>        8888      888E      `888N
`8888L %  ?888   ! u8888cJ888   .d888L .+    .8888Lu=   888E   .u./"888&
 `8888  `-*""   /   "*888*P"    ^"8888*"     ^%888*     888&  d888" Y888*"
   "888.      :"      "Y"          "Y"         "Y"      R888" ` "Y   Y"
     `""***~"`                                           ""
                             https://cortix.org                              
_____________________________________________________________________________
In [3]:
# Create the application network

# Vortex module (single)
vortex = Vortex()
swirl_net.module(vortex)
vortex.show_time = (True,1*const.minute)
vortex.end_time = end_time
vortex.time_step = time_step

for i in range(n_droplets):
    
    # Droplet modules (multiple)
    droplet = Droplet()
    swirl_net.module(droplet)
    droplet.end_time = end_time
    droplet.time_step = time_step
    droplet.bounce = False    # allow droplets to bounce off the ground
    droplet.slip = False          # allow droplets to slip on the ground (otherwise will stick)
    droplet.save = True         # Save the simulation data for post-processing
    
    swirl_net.connect( [droplet, 'external-flow'], [vortex, vortex.get_port('fluid-flow:{}'.format(i))], 'bidirectional')

Verify the network connectivity

In [4]:
# View the Cortix network created

swirl_net.draw()
Out[4]:
network-0 0 Vortex 1 Droplet 0->1 2 Droplet 0->2 3 Droplet 0->3 4 Droplet 0->4 5 Droplet 0->5 6 Droplet 0->6 7 Droplet 0->7 8 Droplet 0->8 9 Droplet 0->9 10 Droplet 0->10 11 Droplet 0->11 12 Droplet 0->12 13 Droplet 0->13 14 Droplet 0->14 15 Droplet 0->15 1->0 2->0 3->0 4->0 5->0 6->0 7->0 8->0 9->0 10->0 11->0 12->0 13->0 14->0 15->0

Run network simulation

In [5]:
# Run the simulation!

swirl.run()
[65266] 2019-08-28 09:54:47,103 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.vortex.Vortex object at 0x152178ba58>
[65266] 2019-08-28 09:54:47,111 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x152178ba20>
[65266] 2019-08-28 09:54:47,121 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217b30f0>
[65266] 2019-08-28 09:54:47,114 - cortix - INFO - Vortex::time[min] = 0.0
[65266] 2019-08-28 09:54:47,132 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217c6908>
[65266] 2019-08-28 09:54:47,144 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217c6ef0>
[65266] 2019-08-28 09:54:47,155 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217cf358>
[65266] 2019-08-28 09:54:47,166 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217cf940>
[65266] 2019-08-28 09:54:47,178 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217cff28>
[65266] 2019-08-28 09:54:47,192 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217d6550>
[65266] 2019-08-28 09:54:47,203 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217d6b38>
[65266] 2019-08-28 09:54:47,215 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217dd128>
[65266] 2019-08-28 09:54:47,227 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217c6cc0>
[65266] 2019-08-28 09:54:47,238 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217dd748>
[65266] 2019-08-28 09:54:47,250 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217ddd30>
[65266] 2019-08-28 09:54:47,262 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217e4358>
[65266] 2019-08-28 09:54:47,272 - cortix - INFO - Launching Module <cortix.examples.droplet_swirl.droplet.Droplet object at 0x15217e4940>
[65266] 2019-08-28 09:54:56,430 - cortix - INFO - Vortex::time[min] = 1.0
[65266] 2019-08-28 09:55:05,939 - cortix - INFO - Vortex::time[min] = 2.0
[65266] 2019-08-28 09:55:14,861 - cortix - INFO - run()::Elapsed wall clock time [s]: 40.39

Results inspection through Cortix

All droplets

In [12]:
'''All droplets' trajectory'''

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

position_histories = list()
for m in swirl_net.modules[1:]:
    position_histories.append( m.liquid_phase.get_quantity_history('position')[0].value )
    
fig = plt.figure(1)
ax = fig.add_subplot(111,projection='3d')
ax.set_title('Droplet Trajectories')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
for p in position_histories:
    x = [u[0] for u in p]
    y = [u[1] for u in p]
    z = [u[2] for u in p]
    ax.plot(x,y,z)
plt.rcParams['figure.figsize'] = [10,10]
plt.show()

Trajectories of all droplets released from random positions at 500-m altitude. Multiprocessing parallel run with number of processes corresponding to all modules: all Droplet modules, 1 Vortex module, and 1 Cortix master process.

In [14]:
'''All droplets' speed'''

fig = plt.figure(2)
plt.xlabel('Time [min]')
plt.ylabel('Speed [m/s]')
plt.title('All Droplets')

for m in swirl_net.modules[1:]:
    speed = m.liquid_phase.get_quantity_history('speed')[0].value
    plt.plot(list(speed.index/60), speed.tolist())

plt.rcParams['figure.figsize'] = [8,5]
plt.grid()
plt.show()
In [8]:
'''All droplets' radial position'''

fig = plt.figure(3)
plt.xlabel('Time [min]')
plt.ylabel('Radial Position [m]')
plt.title('All Droplets')

for m in swirl_net.modules[1:]:
    radial_pos = m.liquid_phase.get_quantity_history('radial-position')[0].value
    plt.plot(list(radial_pos.index/60), radial_pos.tolist())

plt.rcParams['figure.figsize'] = [8,5]
plt.grid()
plt.show()

Individual droplet modules

In [9]:
'''Droplet 1 Module State'''

droplet_1 = swirl_net.modules[2]
(speed_quant,time_unit) = droplet_1.liquid_phase.get_quantity_history('speed')
print('time unit = ',time_unit)
print('speed unit = ',speed_quant.unit)
speed = speed_quant.value
speed.plot(title='Droplet 1: Speed vs Time')
plt.rcParams['figure.figsize'] = [8,5]
plt.grid()
plt.show()
time unit =  s
speed unit =  m/s
In [10]:
(radial_position_quant,time_unit) = droplet_1.liquid_phase.get_quantity_history('radial-position')
print('time unit = ',time_unit)
print('rad. pos. unit = ',radial_position_quant.unit)
rad_pos = radial_position_quant.value
rad_pos.plot(title='Droplet 1: Radial Position vs Time')
plt.rcParams['figure.figsize'] = [8,5]
plt.grid()
plt.show()
time unit =  s
rad. pos. unit =  m

Droplet at Scale

Below is a graph of the number of droplets in the system vs the elapsed time of the corresponding simulation. This scaling experiment uses randomly positioned droplets at altitude of 500 m. Therefore the work for each run is not exactly the same since time integration of the trajectories use time adaptivity and some trajectories take more work to complete than others. However the plot below shows the trend of the one-to-all communication bottlenect with increasing number of droplet modules.

Data collected courtesy of the Idaho National Labs HPC center (https://hpc.inl.gov)

In [11]:
import matplotlib.pyplot as plt

# Droplet number vs Simulation runtimes 
droplet_run_times = [(250, 127), (500, 168), (1000, 346), (2000, 1660), (3000, 2659.24)]

# Calculate the average runtime per droplet
drops_per_sec = [x/y for (x,y) in droplet_run_times]
avg_drop_per_sec = sum(drops_per_sec) / len(drops_per_sec)
print("Average number of droplets handled per second: %.2f" % avg_drop_per_sec)

plt.plot([x for (x, y) in droplet_run_times], [y for (x, y) in droplet_run_times])
plt.title("Wall-Clock Time vs Number of Droplets")
plt.xlabel("Number of Droplets (MPI processes or Cores)")
plt.ylabel("Wall-Clock Time (s)")
plt.grid()
plt.savefig("droplet_trend.png")
plt.show()
Average number of droplets handled per second: 2.03
In [ ]: