#!/usr/bin/env python
# coding: utf-8

# <div style='background-image: url("../../share/images/header.svg") ; padding: 0px ; background-size: cover ; border-radius: 5px ; height: 250px'>
#     <div style="float: right ; margin: 50px ; padding: 20px ; background: rgba(255 , 255 , 255 , 0.7) ; width: 50% ; height: 150px">
#         <div style="position: relative ; top: 50% ; transform: translatey(-50%)">
#             <div style="font-size: xx-large ; font-weight: 900 ; color: rgba(0 , 0 , 0 , 0.8) ; line-height: 100%">Salvus</div>
#             <div style="font-size: large ; padding-top: 20px ; color: rgba(0 , 0 , 0 , 0.5)">2D Global Wave Propagation</div>
#         </div>
#     </div>
# </div>

# In[1]:


import os
import salvus_seismo

# Paths are grabbed from environment variables.
PARAVIEW_BIN = os.environ["PARAVIEW_BIN"]
SALVUS_BIN = os.environ["SALVUS_BIN"]


# In this tutorial we'll run a sesimic simulation on a 2D global Earth model. We'll also use the `salvus_seismo` subpackage to make adding source and receiers a breeze. While we're using a 2D example here so it can be run locally, the steps are the same for running 3D global simulations. First, let's generate a mesh for an Earth model suitable for a 200 second simulation. Use `salvus_mesher`'s command line interface to accomplish this.

# In[2]:


get_ipython().system('python -m pymesher.interface Circular2D --basic.period=200   --chunk2D.max_colatitude=180 --basic.model prem_iso_one_crust --overwrite')


# Now, let's visualize the mesh in Paraview and inspect it.

# In[3]:


get_ipython().system('$PARAVIEW_BIN ./Circular2D_prem_iso_one_crust_200.e')


# Below, we'll rely on `salvus_seismo` to generate a call to Salvus for us. This is a `python` package which makes it simple to generate sources and receivers for use in global 2D/3D Earth models. No more messing around on the command line line as before! See an example below for putting an explosive source near the surface. Check out the comments below for an explanation of the individual parameters.

# In[4]:


import salvus_seismo

# Add a source 1 meter under the surface. Specialize a 2-D "moment tensor", with
# a dominant period of 250 seconds.
src = salvus_seismo.Source(longitude=0.0, depth_in_m=1,
                           m_rr=1e20, m_rp=0, m_pp=1e20,
                           center_frequency=0.004)

# Place a single receiver at a lonitude of 45 degrees. In the 2-D globe, the
# latitude component is zero. Give it a made-up station name.
rec = salvus_seismo.Receiver(longitude=45, depth_in_m=0,
                             network="XX", station="AA")

# Generate the configuration object for salvus_seismo. Note the presence of "salvus_call"
# here. This can be use to pass any special commands, such as mpirun.
# This is a short run for speed reasons - choose a larger endtime for nicer
# visualizations!
config = salvus_seismo.Config(
    mesh_file="./Circular2D_prem_iso_one_crust_200.e",
    end_time=200,
    salvus_call=SALVUS_BIN,
    movie_file_name="test_movie.h5",
    movie_fields=["u_ELASTIC"],
    polynomial_order=4,
    dimensions=2,
    verbose=False)

# Ensure a clean directory. Salvus_seismo will fail to produce the configuration files if
# the directory already exists.
get_ipython().system('rm -rf basic_example/')

# Automatically generate the command line call for salvus. This takes the mesh file as an
# argument, as salvus_seismo will ensure that the source / receiver positions are adjusted
# to respect the mesh topology.
salvus_seismo.generate_cli_call(
    source=src, receivers=[rec], config=config,
    output_folder="basic_example",
    exodus_file="./Circular2D_prem_iso_one_crust_200.e")


# This call should generate a directory with some parameter files inside. All the receiver information is now contained within the file `receivers.toml`, and a snippet containing the command line parameters can be seen in `run_salvus.sh`. We'll take a look at them below.

# In[5]:


get_ipython().system('echo "\\n\\nreceivers.toml:\\n"')
get_ipython().system('cat basic_example/receivers.toml')

get_ipython().system('echo "\\n\\nSalvus run command:\\n"')
get_ipython().system('cat basic_example/run_salvus.sh')


# This run command should contain everything you need to run `Salvus`. You can go ahead and run it as below. Unfortunately, the progress reports in the notebook doesn't seem to work properly (the output is not captured until after the code is finishsed running). You can see the output if you run `sh ./basic_examples/run_salvus` from the command line. Either way, this will run the 2D global simulation.

# In[6]:


get_ipython().system('sh ./basic_example/run_salvus.sh')


# When the simulation is done, we can go ahead and view the movie using the topology output by petsc.

# In[7]:


get_ipython().system('python ./petsc_gen_xdmf.py ./test_movie.h5')


# In[8]:


get_ipython().system('$PARAVIEW_BIN ./test_movie.xmf')


# You don't have to stop here. The source and receivers can parse any ObsPy objects and thus enable a data-driven generation of the Salvus input files. Aside from being convenient this also avoid common pitfalls like wrong unit conversions and manual human intervention.
# 
# This example directly acquires stations coordinates from the IRIS data center. Note that in 2D the latitude will be ignored and the simulation thus happens on the equatorial plane. Moment tensors and vectorial sources will consequently ignore the t (theta, south) components.
# 
# In this case the `generate_cli_call()` takes an additional `exodus_file` argument that takes a mesh. If given, the receivers will be placed exactly relative to the actual mesh surface.
# 
# Also note that usage in 3D is completely identical.

# In[9]:


from obspy.clients.fdsn import Client

src = salvus_seismo.Source(longitude=0.0, depth_in_m=1,
                           m_rr=1e20, m_rp=0, m_pp=1e20,
                           center_frequency=0.004)

c = Client("IRIS")
recs = salvus_seismo.Receiver.parse(c.get_stations(network="IU"))

# Again a short simulation for speed reasons.
config = salvus_seismo.Config(
    mesh_file="./Circular2D_prem_iso_one_crust_200.e",
    end_time=200,
    salvus_call=SALVUS_BIN,
    polynomial_order=4,
    verbose=False,
    dimensions=2)

# Ensure a clean directory.
get_ipython().system('rm -rf webservice_example/')
salvus_seismo.generate_cli_call(
    source=src, receivers=recs, config=config,
    output_folder="webservice_example",
    exodus_file="./Circular2D_prem_iso_one_crust_200.e")


# In[10]:


get_ipython().system('echo "\\n\\nreceivers.toml:\\n"')
get_ipython().system('head -n 10 webservice_example/receivers.toml')

get_ipython().system('echo "\\n\\nSalvus run command:\\n"')
get_ipython().system('cat webservice_example/run_salvus.sh')


# In[11]:


# You can uncomment this line to visualize the receivers on the mesh in Paraview.
# !$PARAVIEW_BIN ./webservice_example/receivers_paraview.csv
get_ipython().system('sh ./webservice_example/run_salvus.sh')


# This has now generated an `hdf5` file which can be opened and visualized with the [ASDF sextant](https://github.com/SeismicData/asdf_sextant). Use it to open `receivers.h5`.