.. _tutorial:

========================
A more detailed tutorial
========================

In the previous tutorial (:doc:`circular_patch_simple`) we provided a high
level overview of the PySPH framework.  No details were provided on equations,
integrators and solvers.  This tutorial assumes that you have read the
previous one.

Recall that in the previous tutorial, a circular patch of fluid with a given
initial velocity field was simulated using a weaky-compressible SPH scheme.
In that example, a ``WCSPHScheme`` object was created in the ``create_scheme``
method.  The details of what exactly the scheme does was not discussed.  This
tutorial explains some of those details by solving the same problem using a
lower-level approach where the actual SPH equations, the integrator, and the
solver are created manually.  This should help a user write their own schemes
or modify an existing scheme.  The full code for this example can be seen in `elliptical_drop_no_scheme.py
<https://github.com/pypr/pysph/tree/master/pysph/examples/elliptical_drop_no_scheme.py>`_.


Imports
~~~~~~~~~~~~~

This example requires a few more imports than the previous case.

the first several lines are imports of various modules:

.. code-block:: python

    import os
    from numpy import array, ones_like, mgrid, sqrt

    # PySPH base and carray imports
    from pysph.base.utils import get_particle_array_wcsph
    from pysph.base.kernels import Gaussian

    # PySPH solver and integrator
    from pysph.solver.application import Application
    from pysph.solver.solver import Solver
    from pysph.sph.integrator import EPECIntegrator
    from pysph.sph.integrator_step import WCSPHStep

    # PySPH sph imports
    from pysph.sph.equation import Group
    from pysph.sph.basic_equations import XSPHCorrection, ContinuityEquation
    from pysph.sph.wc.basic import TaitEOS, MomentumEquation


.. note::

    This is common for all examples that do not use a scheme and it is worth
    noting the pattern of the PySPH imports. Fundamental SPH constructs like
    the kernel and particle containers are imported from the ``base``
    subpackage. The framework related objects like the solver and integrator
    are imported from the ``solver`` subpackage. Finally, we import from the
    ``sph`` subpackage, the physics related part for this problem.

The methods defined for creating the particles are the same as in the previous
tutorial with the exception of the call to
``self.scheme.setup_properties([pa])``.  In this example, we do not create a
scheme, we instead create all the required PySPH objects from the
application.  We do not override the ``create_scheme`` method but instead have
two other methods called ``create_solver`` and ``create_equations`` which
handle this.


Setting up the PySPH framework
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

As we move on, we encounter instantiations of the PySPH framework objects.
These are the :py:class:`pysph.solver.application.Application`,
:py:class:`pysph.sph.integrator.TVDRK3Integrator` and
:py:class:`pysph.solver.solver.Solver` objects.  The ``create_solver`` method
constructs a ``Solver`` instance and returns it as seen below:

.. code-block:: python

        def create_solver(self):
            kernel = Gaussian(dim=2)

            integrator = EPECIntegrator( fluid=WCSPHStep() )

            dt = 5e-6; tf = 0.0076
            solver = Solver(kernel=kernel, dim=2, integrator=integrator,
                            dt=dt, tf=tf, adaptive_timestep=True,
                            cfl=0.05, n_damp=50,
                            output_at_times=[0.0008, 0.0038])

            return solver

As can be seen, various options are configured for the solver, including
initial damping etc.

.. py:currentmodule:: pysph.sph.integrator

Intuitively, in an SPH simulation, the role of the
:py:class:`EPECIntegrator` should be obvious. In the code, we see that we
ask for the "fluid" to be stepped using a :py:class:`WCSPHStep` object. Taking
a look at the ``create_particles`` method once more, we notice that the
**ParticleArray** representing the circular patch was named as `fluid`. So
we're essentially asking the PySPH framework to step or *integrate* the
properties of the **ParticleArray** fluid using :py:class:`WCSPHStep`. It is
safe to assume that the framework takes the responsibility to call this
integrator at the appropriate time during a time-step.

.. py:currentmodule:: pysph.solver.solver

The :py:class:`Solver` is the main driver for the problem. It marshals a
simulation and takes the responsibility (through appropriate calls to the
integrator) to update the solution to the next time step. It also handles
input/output and computing global quantities (such as minimum time step) in
parallel.

Specifying the interactions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

At this stage, we have the particles (represented by the fluid
**ParticleArray**) and the framework to integrate the solution and
marshall the simulation. What remains is to define how to actually go
about updating properties *within* a time step. That is, for each
particle we must "do something". This is where the *physics* for the
particular problem comes in.

For SPH, this would be the pairwise interactions between particles. In PySPH,
we provide a specific way to define the sequence of interactions which is a
*list* of **Equation** objects (see :doc:`../reference/equations`). For the
circular patch test, the sequence of interactions is relatively
straightforward:

    - Compute pressure from the Equation of State (EOS):  :math:`p = f(\rho)`
    - Compute the rate of change of density: :math:`\frac{d\rho}{dt}`
    - Compute the rate of change of velocity (accelerations): :math:`\frac{d\boldsymbol{v}}{dt}`
    - Compute corrections for the velocity (XSPH): :math:`\frac{d\boldsymbol{x}}{dt}`

Care must be taken that the EOS equation should be evaluated for all the
particles before the other equations are evaluated.

.. py:currentmodule:: pysph.sph.equation


We request this in PySPH by creating a list of :py:class:`Equation` instances
in the ``create_equations`` method:

.. code-block:: python

    def create_equations(self):
        equations = [
            Group(equations=[
                TaitEOS(dest='fluid', sources=None, rho0=self.ro,
                        c0=self.co, gamma=7.0),
            ], real=False),

            Group(equations=[
                ContinuityEquation(dest='fluid',  sources=['fluid',]),

                MomentumEquation(dest='fluid', sources=['fluid'],
                                 alpha=self.alpha, beta=0.0, c0=self.co),

                XSPHCorrection(dest='fluid', sources=['fluid']),

            ]),
        ]
        return equations

Each ``Group`` instance is completed before the next is taken up.  Each group
contains a list of ``Equation`` objects.  Each *interaction* is specified
through an :py:class:`Equation` object, which is instantiated with the general
syntax:

.. code-block:: python

   Equation(dest='array_name', sources, **kwargs)

The ``dest`` argument specifies the *target* or *destination*
**ParticleArray** on which this interaction is going to operate
on. Similarly, the ``sources`` argument specifies a *list* of
**ParticleArrays** from which the contributions are sought. For some
equations like the EOS, it doesn't make sense to define a list of
sources and a ``None`` suffices. The specification basically tells PySPH
that for one time step of the calculation:

    - Use the Tait's EOS to update the properties of the fluid array
    - Compute :math:`\frac{d\rho}{dt}` for the fluid from the fluid
    - Compute accelerations for the fluid from the fluid
    - Compute the XSPH corrections for the fluid, using fluid as the source

.. note::

   Notice the use of the **ParticleArray** name "fluid". It is the
   responsibility of the user to ensure that the equation
   specification is done in a manner consistent with the creation of
   the particles.

With the list of equations, our problem is completely defined. PySPH
now knows what to do with the particles within a time step. More
importantly, this information is enough to generate code to carry out
a complete SPH simulation.  For more details on how new equations can be
written please read :ref:`design_overview`.

The example may be run the same way as the previous example::

  $ pysph run elliptical_drop_no_scheme

The resulting output can be analyzed or viewed the same way as in the previous
example.

In the previous example (:doc:`circular_patch_simple`), the equations and solver
are created automatically by the ``WCSPHScheme``. If the ``create_scheme`` is
overwritten and returns a scheme, the ``create_equations`` and ``create_solver``
need not be implemented. For more details on the various application methods,
please see :py:class:`pysph.solver.application.Application`. Implementing other
schemes can be done by either implementing the equations directly as done in
this example or one could implement a new :py:class:`pysph.sph.scheme.Scheme`.