# A more detailed tutorial

In the previous tutorial (Learning the ropes) 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.

## Imports

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

the first several lines are imports of various modules:

    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 `pysph.solver.application.Application`, `pysph.sph.integrator.TVDRK3Integrator` and `pysph.solver.solver.Solver` objects. The `create_solver` method constructs a `Solver` instance and returns it as seen below:

    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.

Intuitively, in an SPH simulation, the role of the `EPECIntegrator` should be obvious. In the code, we see that we ask for the “fluid” to be stepped using a `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 `WCSPHStep`. It is safe to assume that the framework takes the responsibility to call this integrator at the appropriate time during a time-step.

The `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 SPH equations). For the circular patch test, the sequence of interactions is relatively straightforward:

> - Compute pressure from the Equation of State (EOS): p = f(\rho)
> - Compute the rate of change of density: \frac{d\rho}{dt}
> - Compute the rate of change of velocity (accelerations): \frac{d\boldsymbol{v}}{dt}
> - Compute corrections for the velocity (XSPH): \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.

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

    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 `Equation` object, which is instantiated with the general syntax:

    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 \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 The PySPH framework.

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 (Learning the ropes), 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. Implementing other schemes can be done by either implementing the equations directly as done in this example or one could implement a new `Scheme`.