"""Dam break past an obstacle using initial data generated by SPHysics

For benchmarking, we use the input geometry and discretization as the
SPHYSICS Case 5
(https://wiki.manchester.ac.uk/sphysics/index.php/SPHYSICS_Home_Page)

We only require the INDAT and IPART files generated by SPHysics. These
define respectively, the numerical parameters and the initial particle
data used for the run. The rest of the problem is set-up in the usual
way.

"""
import numpy

from pysph.sph.equation import Group
from pysph.base.kernels import CubicSpline, Gaussian
from pysph.sph.wc.basic import TaitEOS, TaitEOSHGCorrection, MomentumEquation
from pysph.sph.basic_equations import ContinuityEquation, XSPHCorrection

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

from pysph.tools.sphysics import sphysics2pysph

# problem dimensionality
dim = 3

# suggested initial time step and final time
dt = 1e-5
tf = 2.0

# physical constants for the run loaded from SPHysics INDAT
indat = numpy.loadtxt('INDAT.gz')
H = float( indat[10] )
B = float( indat[11] )
gamma = float( indat[12] )
eps = float( indat[14] )
rho0 = float( indat[15] )
alpha = float( indat[16] )
beta = 0.0
c0 = numpy.sqrt( B*gamma/rho0 )

# particle factory function
def create_particles(ipart='IPART.gz', indat='INDAT.gz', **kwargs):
     return sphysics2pysph(ipart, indat, vtk=False)

# Create the application.
app = Application()

# Specify the kernel
kernel = CubicSpline(dim=3)

# Setup the integrator.
integrator = EPECIntegrator(fluid=WCSPHStep(),boundary=WCSPHStep())

# Create a solver. damping time is taken as 0.1% of the final time
solver = Solver(dim=dim, kernel=kernel, integrator=integrator,
                adaptive_timestep=True, tf=tf, dt=dt, n_damp=100)

# create the equations
equations = [

    # Equation of state
    Group(equations=[

            TaitEOS(dest='fluid', sources=None, rho0=rho0, c0=c0, gamma=gamma),
            TaitEOSHGCorrection(dest='boundary', sources=None, rho0=rho0, c0=c0, gamma=gamma),

            ], real=False),

    # Continuity Momentum and XSPH equations
    Group(equations=[
               ContinuityEquation(dest='fluid', sources=['fluid', 'boundary']),
               ContinuityEquation(dest='boundary', sources=['fluid']),

               MomentumEquation(dest='fluid', sources=['fluid', 'boundary'],
                                alpha=alpha, beta=beta, gz=-9.81,
                                tensile_correction=True),

               # Position step with XSPH
               XSPHCorrection(dest='fluid', sources=['fluid'], eps=eps)
               ])
    ]


# Setup the application and solver.  This also generates the particles.
app.setup(solver=solver, equations=equations,
          particle_factory=create_particles)

app.run()
