#!/usr/bin/env python2.7
# -*- coding: UTF-8 -*-
#
# Copyright (c) 2014, Yung-Yu Chen <yyc@solvcon.net>
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
# - Neither the name of the SOLVCON nor the names of its contributors may be
#   used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.


"""
This is an example for solving the problem of oblique shock reflection.
"""


import os # Python standard library
import numpy as np # http://www.numpy.org
import solvcon as sc # SOLVCON
from solvcon.parcel import gas # A specific SOLVCON solver package we'll use


class ObliqueShockReflection(gas.ObliqueShockRelation):
    def __init__(self, gamma, theta, mach1, rho1, p1, T1):
        super(ObliqueShockReflection, self).__init__(gamma=gamma)
        # Angles and Mach numbers.
        self.theta = theta
        self.mach1 = mach1
        self.beta1 = beta1 = self.calc_shock_angle(mach1, theta)
        self.mach2 = mach2 = self.calc_dmach(mach1, beta1)
        self.beta2 = beta2 = self.calc_shock_angle(mach2, theta)
        self.mach3 = mach3 = self.calc_dmach(mach2, beta2)
        # Flow properties in the first zone.
        self.rho1 = rho1
        self.p1 = p1
        self.T1 = T1
        self.a1 = np.sqrt(gamma*p1/rho1)
        # Flow properties in the second zone.
        self.rho2 = rho2 = rho1 * self.calc_density_ratio(mach1, beta1)
        self.p2 = p2 = p1 * self.calc_pressure_ratio(mach1, beta1)
        self.T2 = T2 = T1 * self.calc_temperature_ratio(mach1, beta1)
        self.a2 = np.sqrt(gamma*p2/rho2)
        # Flow properties in the third zone.
        self.rho3 = rho3 = rho2 * self.calc_density_ratio(mach2, beta2)
        self.p3 = p3 = p2 * self.calc_pressure_ratio(mach2, beta2)
        self.T3 = T3 = T2 * self.calc_temperature_ratio(mach2, beta2)
        self.a3 = np.sqrt(gamma*p3/rho3)

    def __str__(self):
        msg = 'Relation of reflected oblique shock:\n'
        msg += '- theta = %5.2f deg (flow angle)\n' % (self.theta/np.pi*180)
        msg += '- beta1 = %5.2f deg (shock angle)\n' % (self.beta1/np.pi*180)
        msg += '- beta1 = %5.2f deg (shock angle)\n' % (self.beta2/np.pi*180)
        def property_string(zone):
            values = [getattr(self, '%s%d' % (key, zone))
                      for key in 'mach', 'rho', 'p', 'T', 'a']
            messages = [' %6.3f' % val for val in values]
            return ''.join(messages)
        msg += '- mach, rho, p, T, a (1) =' + property_string(1) + '\n'
        msg += '- mach, rho, p, T, a (2) =' + property_string(2) + '\n'
        msg += '- mach, rho, p, T, a (3) =' + property_string(3)
        return msg

    @property
    def hookcls(self):
        relation = self
        class _ShowRelation(sc.MeshHook):
            def preloop(self):
                for msg in str(relation).split('\n'):
                    self.info(msg + '\n')
            postloop = preloop
        return _ShowRelation


class RectangleMesher(object):
    """
    Representation of a rectangle and the Gmsh meshing helper.

    :ivar lowerleft: (x0, y0) of the rectangle.
    :type lowerleft: tuple
    :ivar upperright: (x1, y1) of the rectangle.
    :type upperright: tuple
    :ivar edgelength: Length of each mesh cell edge.
    :type edgelength: float
    """

    GMSH_SCRIPT = """
    // vertices.
    Point(1) = {%(x1)g,%(y1)g,0,%(edgelength)g};
    Point(2) = {%(x0)g,%(y1)g,0,%(edgelength)g};
    Point(3) = {%(x0)g,%(y0)g,0,%(edgelength)g};
    Point(4) = {%(x1)g,%(y0)g,0,%(edgelength)g};
    // lines.
    Line(1) = {1,2};
    Line(2) = {2,3};
    Line(3) = {3,4};
    Line(4) = {4,1};
    // surface.
    Line Loop(1) = {1,2,3,4};
    Plane Surface(1) = {1};
    // physics.
    Physical Line("upper") = {1};
    Physical Line("left") = {2};
    Physical Line("lower") = {3};
    Physical Line("right") = {4};
    Physical Surface("domain") = {1};
    """.strip()

    def __init__(self, lowerleft, upperright, edgelength):
        assert 2 == len(lowerleft)
        self.lowerleft = tuple(float(val) for val in lowerleft)
        assert 2 == len(upperright)
        self.upperright = tuple(float(val) for val in upperright)
        self.edgelength = float(edgelength)

    def __call__(self, cse):
        x0, y0 = self.lowerleft
        x1, y1 = self.upperright
        script_arguments = dict(
            edgelength=self.edgelength, x0=x0, y0=y0, x1=x1, y1=y1)
        cmds = self.GMSH_SCRIPT % script_arguments
        cmds = [cmd.rstrip() for cmd in cmds.strip().split('\n')]
        gmh = sc.Gmsh(cmds)()
        return gmh.toblock(bcname_mapper=cse.condition.bcmap)


def generate_bcmap(relation):
    # Flow properties (fp).
    fpb = {
        'gamma': relation.gamma, 'rho': relation.rho1,
        'v2': 0.0, 'v3': 0.0, 'p': relation.p1,
    }
    fpb['v1'] = relation.mach1*np.sqrt(relation.gamma*fpb['p']/fpb['rho'])
    fpt = fpb.copy()
    fpt['rho'] = relation.rho2
    fpt['p'] = relation.p2
    V2 = relation.mach2 * relation.a2
    fpt['v1'] = V2 * np.cos(relation.theta)
    fpt['v2'] = -V2 * np.sin(relation.theta)
    fpi = fpb.copy()
    # BC map.
    bcmap = {
        'upper': (sc.bctregy.GasInlet, fpt,),
        'left': (sc.bctregy.GasInlet, fpb,),
        'right': (sc.bctregy.GasNonrefl, {},),
        'lower': (sc.bctregy.GasWall, {},),
    }
    return bcmap


class ReflectionProbe(gas.ProbeHook):
    """
    Place a probe for the flow properties in the reflected region.
    """

    def __init__(self, cse, **kw):
        """
        :param relation: Instruct shock relations.
        :type relation: ObliqueShockReflection
        :param rect: Instruct the mesh rectangle.
        :type rect: RectangleMesher
        """
        rect = kw.pop('rect')
        self.relation = relation = kw.pop('relation')
        factor = kw.pop('factor', 0.9)
        # Detemine location.
        theta = relation.theta
        beta1 = relation.beta1
        beta2 = relation.beta2
        x0, y0 = rect.lowerleft
        x1, y1 = rect.upperright
        lgh = (y1-y0) / np.tan(beta1)
        hgt = factor * (x1-x0-lgh) * np.tan((beta2-theta)/2)
        lgh = factor * (x1-x0-lgh) + lgh
        poi = ('poi', x0+lgh, y0+hgt, 0.0)
        # Call super.
        kw['coords'] = (poi,)
        kw['speclst'] = ['M', 'rho', 'p']
        super(ReflectionProbe, self).__init__(cse, **kw)

    def postloop(self):
        super(ReflectionProbe, self).postloop()
        rel = self.relation
        self.info('Probe result at %s:\n' % self.points[0])
        M, rho, p = self.points[0].vals[-1][1:]
        self.info('- mach3 = %.3f/%.3f (error=%%%.2f)\n' % (
            M, rel.mach3, abs((M-rel.mach3)/rel.mach3)*100))
        self.info('- rho3  = %.3f/%.3f (error=%%%.2f)\n' % (
            rho, rel.rho3, abs((rho-rel.rho3)/rel.rho3)*100))
        self.info('- p3    = %.3f/%.3f (error=%%%.2f)\n' % (
            p, rel.p3, abs((p-rel.p3)/rel.p3)*100))


def obrf_base(
    casename=None, psteps=None, ssteps=None, edgelength=None,
    gamma=None, density=None, pressure=None, mach=None, theta=None, **kw):
    """
    Base configuration of the simulation and return the case object.

    :return: The created Case object.
    :rtype: solvcon.parcel.gas.GasCase
    """

    ############################################################################
    # Step 1: Obtain the analytical solution.
    ############################################################################

    # Calculate the flow properties in all zones separated by the shock.
    relation = ObliqueShockReflection(gamma=gamma, theta=theta, mach1=mach,
                                      rho1=density, p1=pressure, T1=1)

    ############################################################################
    # Step 2: Instantiate the simulation case.
    ############################################################################

    # Create the mesh generator.  Keep it for later use.
    mesher = RectangleMesher(lowerleft=(0,0), upperright=(4,1),
                             edgelength=edgelength)
    # Set up case.
    cse = gas.GasCase(
        # Mesh generator.
        mesher=mesher,
        # Mapping boundary-condition treatments.
        bcmap=generate_bcmap(relation),
        # Use the case name to be the basename for all generated files.
        basefn=casename,
        # Use `cwd`/result to store all generated files.
        basedir=os.path.abspath(os.path.join(os.getcwd(), 'result')),
        # Debug and capture-all.
        debug=False, **kw)

    ############################################################################
    # Step 3: Set up delayed callbacks.
    ############################################################################

    # Field initialization and derived calculations.
    cse.defer(gas.FillAnchor, mappers={'soln': gas.GasSolver.ALMOST_ZERO,
                                       'dsoln': 0.0, 'amsca': gamma})
    cse.defer(gas.DensityInitAnchor, rho=density)
    cse.defer(gas.PhysicsAnchor, rsteps=ssteps)
    # Report information while calculating.
    cse.defer(relation.hookcls)
    cse.defer(gas.ProgressHook, linewidth=ssteps/psteps, psteps=psteps)
    cse.defer(gas.CflHook, fullstop=False, cflmax=10.0, psteps=ssteps)
    cse.defer(gas.MeshInfoHook, psteps=ssteps)
    cse.defer(ReflectionProbe, rect=mesher, relation=relation, psteps=ssteps)
    # Store data.
    cse.defer(gas.PMarchSave,
              anames=[('soln', False, -4),
                      ('rho', True, 0),
                      ('p', True, 0),
                      ('T', True, 0),
                      ('ke', True, 0),
                      ('M', True, 0),
                      ('sch', True, 0),
                      ('v', True, 0.5)],
              psteps=ssteps)

    ############################################################################
    # Final: Return the configured simulation case.
    ############################################################################
    return cse


@gas.register_arrangement
def obrf(casename, **kw):
    return obrf_base(
        # Required positional argument for the name of the simulation case.
        casename,
        # Arguments to the base configuration.
        ssteps=200, psteps=4, edgelength=0.1,
        gamma=1.4, density=1.0, pressure=1.0, mach=3.0, theta=10.0/180*np.pi,
        # Arguments to GasCase.
        time_increment=7.e-3, steps_run=600, **kw)


@gas.register_arrangement
def obrf_fine(casename, **kw):
    return obrf_base(
        casename,
        ssteps=200, psteps=4, edgelength=0.02,
        gamma=1.4, density=1.0, pressure=1.0, mach=3.0, theta=10.0/180*np.pi,
        time_increment=1.e-3, steps_run=4000, **kw)


if __name__ == '__main__':
    sc.go()

# vim: set ff=unix fenc=utf8 ft=python ai et sw=4 ts=4 tw=79:
