<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Numerical Wave Tank (3D Version)
This notebook provides a numerical model of a wave tank that would look something like this schematically
Step1:  IPython Parallel "Magics"
As a first example of executing code in parallel on our, we use the cell magic %%px to execute some basic MPI code on every engine in our cluster.
Step2:  Load Proteus
Proteus is a Python package consiting of multiple modules. We are going to define a problem and run a simulation interactively, so we pull in the iproteus (interactive proteus) module to set up a basic environment.
Notes
Step3:  Define the tank geometry
Importing the tank module reads a text file for the bathymetry transect and generates a polygonal domain. There are several options for how we can modify it
Step4:  Physics and Numerics
Load the modules the define the equations to be solved and the numerical methods to use.
Step5:  Numerical Solution Object
No we create an the numerical wavetank object and set it up to run in a thread on each engine.
Step10:  Define some functions to help monitor the calculation
Step11:  Run the tank
Step14:  Post-process the numerical solution
<ASSISTANT_TASK:>
Python Code:
from IPython.parallel import Client, error
cluster = Client(profile='mpi')
view = cluster[:]
Explanation: Numerical Wave Tank (3D Version)
This notebook provides a numerical model of a wave tank that would look something like this schematically:
The numerical model is provided by Proteus (http://proteus.usace.army.mil)
Proteus is being developed at ERDC and includes 2D and 3D versions of this numerical model
This notebook will stick to 2D because it's faster and easier to pre- and post-process.
Mathematical Model
To model a tank with breaking and overtopping, we solve the following air/water flow equations in 2D (a vertical slice):
$$
\begin{eqnarray}
\nabla \cdot \left( \epsilon \mathbf v \right) &=& s \
\frac{\partial \left( \epsilon \mathbf v \right)}{\partial t} +\nabla \cdot  \left(\epsilon \mathbf v \otimes \mathbf v\right) -  \nabla \cdot \left{\epsilon \nu_t \left(\nabla  \mathbf v+\nabla \mathbf v^T\right) \right} + \mathbf r + \frac{\epsilon \nabla p}{\rho} &=& 0 \
\frac{\partial \phi}{\partial t} + \mathbf v \cdot \nabla \phi &=& 0 \
\|\nabla \phi\| &=& 1 \
\frac{\partial \theta}{\partial t} + \nabla \cdot \left( \theta \epsilon \mathbf v \right) &=& 0 \
H(\phi) &=& \theta 
\end{eqnarray}
$$
Notes about this demo
This software can be run locally on mac, windows, and linux if the Proteus stack is installed.
It can also be accessed with a cloud service like Sagemath Cloud, Wakari, or an internal cloud server if we invest in setting one up (you are probably viewing on Sagemath Cloud)
This particular notebook can also be saved as a Python script and run on HPC machines if they have the Proteus software stack installed
Grid generation, specification of input parameters like waves, boundary conditions, and material properties can all be viewed, set from pre-existing definitions, or redefined in this notebook.
The simulation is kicked off in the notebook as a seperate thread and monitored
The data is archived  in XDMF, which can be opened upon completion of the simulation for post-processing, furthermore additional simulations and multiple scripted runs can be orchestrated by modifying this notebook.
Parallelism
We need to use multiple processors to compute the solution efficiently in 3D and sometimes even in 2D. 
The IPython infrastructure allows us to attach to a set of "engines" (e.g. MPI tasks on an HPC machine)
End of explanation
%%px
from mpi4py import MPI
mpi = MPI.COMM_WORLD
bcast = mpi.bcast
barrier = mpi.barrier
rank = mpi.rank
print "MPI rank: %i/%i" % (mpi.rank,mpi.size)
Explanation: IPython Parallel "Magics"
As a first example of executing code in parallel on our, we use the cell magic %%px to execute some basic MPI code on every engine in our cluster.
End of explanation
%%px
import sys
from proteus.iproteus import * 
%%px
from proteus import default_n, default_s, default_so
Profiling.logLevel=5
Profiling.verbose=False
Explanation: Load Proteus
Proteus is a Python package consiting of multiple modules. We are going to define a problem and run a simulation interactively, so we pull in the iproteus (interactive proteus) module to set up a basic environment.
Notes:
Proteus runs with logging that records verying amounts of information by setting logLevel from 1 to 11
The log is stored  in a .log file, which can be downloaded
If Profiling.verbose is set  to True then logging will show up in certain output cells of the notebook, which is usually not what you want
End of explanation
#!tar czvf mydir.tar.gz  ls_p.py tank_so.py twp_navier_stokes_n.py redist_n.py   twp_navier_stokes_p.py redist_p.py vof_n.py vof_p.py ls_consrv_n.py ls_consrv_p.py ls_n.py tank.py
#!scp mydir.tar.gz spirit01.afrl.hpc.mil:
#%%px 
#import subprocess
#if rank == 0:
  #  status_gz = subprocess.call("gunzip " +  "mydir.tar.gz", shell=True)
    #status_tar = subprocess.call("tar "+"xvf mydir.tar",shell=True)
#    os.system('tar xzvf mydir.tgz')
%%px
import floating_bar,floating_bar_so
Explanation: Define the tank geometry
Importing the tank module reads a text file for the bathymetry transect and generates a polygonal domain. There are several options for how we can modify it:
add a file upload box to just upload a different bottom bathymetry
manipulate an interactive version of the resulting domain (an interactive version of the plot below)
manipulate vertices directly in domain.vertices and domain.segments
End of explanation
%%px
from proteus import Comm
from petsc4py import PETSc
so = floating_bar_so
so.tnList = so.tnList
pList=[]
nList=[]
so.sList=[]
OptDB = PETSc.Options()
for (p,n) in so.pnList:
    so.sList.append(default_s)
    pList.append(__import__(p))
    nList.append(__import__(n))
    pList[-1].name = p
    nList[-1].multilevelLinearSolver = default_n.KSP_petsc4py
    nList[-1].levelLinearSolver = default_n.KSP_petsc4py
    OptDB.setValue(nList[-1].linear_solver_options_prefix+"ksp_type", "preonly")
    OptDB.setValue(nList[-1].linear_solver_options_prefix+"pc_type", "lu")
    OptDB.setValue(nList[-1].linear_solver_options_prefix+"pc_factor_mat_solver_package","superlu_dist")
opts.save_dof = True
opts.dataDir='.'
opts.probDir='.'
opts.logLevel=7
opts.verbose=True
Explanation: Physics and Numerics
Load the modules the define the equations to be solved and the numerical methods to use.
End of explanation
%%px
ns = NumericalSolution.NS_base(so, pList, nList, so.sList, opts)
%%px --noblock
from threading import Thread
# Create a thread wrapper for the simulation.  The target must be an argument-less
# function so we wrap the call to `calculateSolution` in a simple lambda:
simulation_thread = Thread(target = lambda : ns.calculateSolution('run1'))
Explanation: Numerical Solution Object
No we create an the numerical wavetank object and set it up to run in a thread on each engine.
End of explanation
import numpy as np
import numpy
from pythreejs import *
from IPython.core.display import clear_output
def plot_current_results(in_place=True):
    import numpy as np
    #from mpl_toolkits.mplot3d import Axes3D
    #import matplotlib.pyplot as  plt
    Makes a blocking call to retrieve remote data and displays the solution mesh
    as a contour plot.
    
    Parameters
    ----------
    in_place : bool
        By default it calls clear_output so that new plots replace old ones.  Set
        to False to allow keeping of all previous outputs.
    
    global nn,x,y,u,vertices,triangles,domain,rigid_bar
    # We make a blocking call to load the remote data from the simulation into simple named 
    # variables we can read from the engine namespaces
    #load_simulation_globals()
    view.apply_sync(load_simulation_globals)
    # And now we can use the view to read these variables from all the engines.  Then we
    # concatenate all of them into single arrays for local plotting
    #x = np.concatenate(view['x'])
    #y = np.concatenate(view['y'])
    #z = np.concatenate(view['z'])
    vertices = np.concatenate(view['vertices'])
    shifts = np.cumsum([0]+view['nn'][:-1])
    flat_triangles = np.concatenate([ tri + shift for tri,shift in zip(view['triangles'], shifts) ])
    #flat_triangles=triangles
    # We can now call the matplotlib plotting function we need
    #fig, ax = plt.subplots(subplot_kw=dict(aspect='equal'))
    #fig = plt.figure()
    #ax = fig.gca(projection='3d')
    #print x.shape,y.shape,flat_triangles.shape
    #print flat_triangles.flat[:].max(),flat_triangles.flat[:].min()
    #ax.plot_trisurf(X=x, Y=y, TRIANGLES=flat_triangles)#, Z=z)
    #ax.plot_trisurf(x,y,z)
    widget_surface.geometry = FaceGeometry(vertices=list(vertices.flatten()),
                                           face3=list(flat_triangles.flatten()))
    #from matplotlib import tri
    #mesh = tri.Triangulation(x,y,flat_triangles)
    #help(mesh)
    #ip = tri.LinearTriInterpolator(triangulation, u, trifinder=None)
    #X = np.linspace(0,tank.domain.L[0],40)
    #Z = np.linspace(0,tank.domain.L[1],40)
    #U = np.zeros((40,40),'d')
    #W = np.zeros((40,40),'d')
    #U = ip(X,Z)
    # We clear the notebook output before plotting this if in-place plot updating is requested
    #if in_place:
    #    clear_output()
    #display(fig)
    #return fig
    return 0
def load_simulation_globals():
    Put some variables we need in engine namespace.
    These can then be retrieved by clients for inspection, visualization, etc.
    
    global nn, vertices,x, y, z, triangles,domain,rigid_bar
    isosurface = ns.auxiliaryVariables[ns.modelList[2].name][0]
    rigid_bar = ns.auxiliaryVariables[ns.modelList[0].name][0]
    domain = ns.pList[0].domain
    nodes = isosurface.nodes_array
    triangles = isosurface.elements_array
    x = nodes[:,0]
    y = nodes[:,1]
    z = nodes[:,2]
    vertices = nodes
    nn = len(x)
    
def simulation_alive():
    Return True if the simulation thread is still running on any engine.
    
    #return simulation_thread.is_alive()
    return any(view.apply_sync(lambda : simulation_thread.is_alive()))
def monitor_simulation(refresh=5.0, plots_in_place=True):
    Monitor the simulation progress and call plotting routine.
    Supress KeyboardInterrupt exception if interrupted, ensure that the last 
    figure is always displayed and provide basic timing and simulation status.
    Parameters
    ----------
    refresh : float
      Refresh interval between calls to retrieve and plot data.  The default
      is 5s, adjust depending on the desired refresh rate, but be aware that 
      very short intervals will start having a significant impact.
    plots_in_place : bool
       If true, every new figure replaces the last one, producing a (slow)
       animation effect in the notebook.  If false, all frames are plotted
       in sequence and appended in the output area.
    
    import datetime as dt, time
    
    if not simulation_alive():
        plot_current_results(in_place=plots_in_place)
        plt.close('all')
        print 'Simulation has already finished, no monitoring to do.'
        return
    
    t0 = dt.datetime.now()
    fig = None
    try:
        while simulation_alive():
            plot_current_results(in_place=plots_in_place)
            #plt.close('all') # prevent re-plot of old figures
            time.sleep(refresh) # so we don't hammer the server too fast
    except (KeyboardInterrupt):#, error.TimeoutError):
        msg = 'Monitoring interrupted, simulation is ongoing!'
    else:
        msg = 'Simulation completed!'
    tmon = dt.datetime.now() - t0
    #if plots_in_place and fig is not None:
    #    clear_output()
    #    display(fig)
    print msg
    print 'Monitored for: %s.' % tmon
Explanation: Define some functions to help monitor the calculation
End of explanation
%px simulation_thread.start()
view.apply_sync(load_simulation_globals)
vertices = np.concatenate(view['vertices'])
shifts = np.cumsum([0]+view['nn'][:-1])
flat_triangles = np.concatenate([ tri + shift for tri,shift in zip(view['triangles'], shifts) ])
from pythreejs import *
from IPython.display import display
domain = view['domain'][0]
vertices = vertices#verticesview['vertices'][0]
triangles=flat_triangles#view['triangles'][0]
#rigid_bar = view['rigid_bar'][0]
#L = view['ns'][0].pList[0].domain.L
#x_ll = view['ns'][0].pList[0].domain.x_ll
#b = view['ns'][0].nList[0].auxiliaryVariables[0].bar
#sx,sy,sz = bar.boxsize
sx,sy,sz = (0.5,0.5,0.5)
widget_bar = Mesh(geometry=BoxGeometry(width=sx,height=sy,depth=sz),
                    material=LambertMaterial(color='black'),
                    position=(0.5,0.5,0.5))#b.getPosition())
#widget.quaternion_from_rotation(b.getRotation())
center = (0.5*(domain.x[0]+domain.L[0]),
          0.5*(domain.x[1]+domain.L[1]),
          0.5*(domain.x[2]+domain.L[2]))
widget_surface = Mesh(geometry=FaceGeometry(vertices=list(vertices.flatten()),face3=list(flat_triangles.flatten())),
                      material=LambertMaterial(color='aqua'))
#widget.quaternion_from_rotation(b.getRotation())
center = (0.5*(domain.x[0]+domain.L[0]),
          0.5*(domain.x[1]+domain.L[1]),
          0.5*(domain.x[2]+domain.L[2]))
          
children=[widget_bar,widget_surface]
children.append(Mesh(geometry=BoxGeometry(width=domain.L[0],height=domain.L[1],depth=domain.L[2]), 
                     material=LambertMaterial(color=0xffffff, #color=0xccccff, 
                                              opacity=0.7,
                                              refractionRatio=0.985, 
                                              reflectivity= 0.9,
                                              transparent=True),
                     position=center))
#children.append(Mesh(geometry=PlaneGeometry(width=10,height=10),material=BasicMaterial(color='blue')))
children.append(AmbientLight(color=0x777777))
scene = Scene(children=children)
c = PerspectiveCamera(target=center,
                      position=(center[0],center[1]+2*domain.L[1],center[2]), 
                      up=[0,0,1], 
                      children=[DirectionalLight(color='white', 
                                                 position=(center[0],center[1],center[2]+2*domain.L[2]), 
                                                 intensity=0.5)])
renderer = Renderer(camera=c, scene = scene, controls=[OrbitControls(controlling=c)])
renderer.background='gray'
ar = float(renderer.height)/(renderer.width)
renderer.width = 600
renderer.height = int(ar*renderer.width)
display(renderer)
monitor_simulation(refresh=5.0)
Explanation: Run the tank
End of explanation
%%px
from tables import  openFile
archive = openFile('tank_p%d.h5' % (rank,),'r')
def load_post_simulation_globals(it):
    Put some variables we need in engine namespace.
    These can then be retrieved by clients for inspection, visualization, etc.
    
    global phi
    print it
    phi=archive.getNode("/phi"+`it`)[:];
x = np.concatenate(view['x'])
y = np.concatenate(view['y'])
shifts = numpy.cumsum([0]+view['nn'][:-1])
flat_triangles = np.concatenate([ tri + shift for tri,shift in zip(view['triangles'], shifts) ])
triplot(x,y,flat_triangles)
!rm phi*png
import tank
for it,t in enumerate(view['tank_so.tnList'][0]):
    view.apply_sync(load_post_simulation_globals,it)
    phi = np.concatenate(view['phi'])
    pyplot.clf()
    plt.xlabel(r'z[m]')
    plt.ylabel(r'x[m]')
    colors = ['b','g','r','c','m','y','k','w']
    pylab.xlim(tank.domain.x[0]-0.1*tank.domain.L[0],tank.domain.x[0]+tank.domain.L[0]+0.1*tank.domain.L[0])    
    pyplot.axis('equal')
    for si,s in enumerate(tank.segments):
        pyplot.plot([tank.domain.vertices[s[0]][0],
                     tank.domain.vertices[s[1]][0]],
                    [tank.domain.vertices[s[0]][1],
                     tank.domain.vertices[s[1]][1]],
                    color=colors[tank.domain.segmentFlags[si]-1],
                    linewidth=2,
                    marker='o')
    pyplot.tricontour(x,y,flat_triangles,phi,[0])
    pyplot.title('T=%2.2f' % (t,))
    pyplot.savefig('phi%4.4d.png' % (it,))
!avconv -y -i phi%4d.png -c libx264 -qscale 1 tankPhi.mp4
from IPython.core.display import HTML
data_uri_mp4 = open("tankPhi.mp4", "rb").read().encode("base64").replace("\n", "")
video_tag = <video controls>
<source type ="video/mp4" src="data:video/mp4;base64,{mp4}"/>
Your browser does not support the video tag
</video>.format(mp4=data_uri_mp4)
HTML(data=video_tag)
from IPython.display import FileLink,FileLinks
FileLink('tankPhi.mp4')
Explanation: Post-process the numerical solution
End of explanation
<END_TASK>