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

"""
***********************************************************************************
                           tutorial_che_5.py
                DAE Tools: pyDAE module, www.daetools.com
                Copyright (C) Dragan Nikolic
***********************************************************************************
DAE Tools is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License version 3 as published by the Free Software
Foundation. DAE Tools is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with the
DAE Tools software; if not, see <http://www.gnu.org/licenses/>.
************************************************************************************
"""
__doc__ = """
Similar to the chem. eng. example 4, this example shows a comparison between the analytical
results and the discretised population balance equations results solved using the cell
centered finite volume method employing the high resolution upwind scheme (Van Leer
k-interpolation with k = 1/3) and a range of flux limiters.

This tutorial can be run from the console only.

The problem is from the section 4.1.2 Size-independent growth II of the following article:

- Nikolic D.D., Frawley P.J. Application of the Lagrangian Meshfree Approach to
  Modelling of Batch Crystallisation: Part II – An Efficient Solution of Integrated CFD
  and Population Balance Equations. Preprints 2016, 20161100128.
  `doi:10.20944/preprints201611.0012.v1 <http://dx.doi.org/10.20944/preprints201611.0012.v1>`_

and also from the section 3.2 Size-independent growth of the following article:

- Qamar S., Elsner M.P., Angelov I.A., Warnecke G., Seidel-Morgenstern A. (2006)
  A comparative study of high resolution schemes for solving population balances in crystallization.
  Computers and Chemical Engineering 30(6-7):1119-1131.
  `doi:10.1016/j.compchemeng.2006.02.012 <http://dx.doi.org/10.1016/j.compchemeng.2006.02.012>`_

Again, the growth-only crystallisation process was considered with the constant growth rate
of 0.1μm/s and with the different initial number density function:

.. code-block:: none

   n(L,0):                      0, if        L <= 2.0μm
                             1E10, if  2μm < L <= 10μm (region I)
                                0, if 10μm < L <= 18μm
         1E10*cos^2(pi*(L-26)/64), if 18μm < L <= 34μm (region II)
                                0, if 34μm < L <= 42μm
         1E10*sqrt(1-(L-50)^2/64), if 42μm < L <= 58μm (region III)
                                0, if 58μm < L <= 66μm
   1E10*exp(-(L-70)^2/(2sigma^2)), if 66μm < L <= 74μm (region IV)
                                0, if 74μm < L

The crystal size in the range of [0, 100]μm was discretised into 200 elements.
The analytical solution in this case is equal to the initial profile translated right
in time by a distance Gt (the growth rate multiplied by the time elapsed in the process).

Comparison of L1- and L2-norms (ni_HR - ni_analytical):

.. code-block:: none

   -------------------------------------
           Scheme  L1         L2
   -------------------------------------
         superbee  4.464e+10  1.015e+10
            smart  4.727e+10  1.120e+10
            Koren  4.861e+10  1.141e+10
            Sweby  5.435e+10  1.142e+10
               MC  5.129e+10  1.162e+10
           HQUICK  5.531e+10  1.194e+10
             HCUS  5.528e+10  1.194e+10
    vanLeerMinmod  5.600e+10  1.202e+10
          vanLeer  5.814e+10  1.225e+10
            ospre  6.131e+10  1.252e+10
            UMIST  6.181e+10  1.259e+10
            Osher  6.690e+10  1.275e+10
       vanAlbada1  6.600e+10  1.281e+10
           minmod  7.751e+10  1.360e+10
       vanAlbada2  7.901e+10  1.413e+10
   -------------------------------------

The comparison of number density functions between the analytical solution and the
solution obtained using high-resolution scheme with the Superbee flux limiter at t=100s:

.. image:: _static/tutorial_che_5-results.png
   :width: 500px

The comparison of number density functions between the analytical solution and the
solution obtained using high-resolution scheme with the Koren flux limiter at t=100s:

.. image:: _static/tutorial_che_5-results2.png
   :width: 500px
"""

import sys, numpy
from daetools.pyDAE import *
from daetools.pyDAE.data_reporters import *
from time import localtime, strftime
try:
    from .fl_analytical import run_analytical
except:
    from fl_analytical import run_analytical

# Standard variable types are defined in variable_types.py
from pyUnits import m, g, kg, s, K, mol, kmol, J, um

pbm_number_density_t = daeVariableType("pbm_number_density_t", m**(-1), 0.0, 1.0e+40,  0.0, 1e-0)

class modelMoC(daeModel):
    def __init__(self, Name, G, Phi, Parent = None, Description = ""):
        daeModel.__init__(self, Name, Parent, Description)

        self.G   = Constant(G)
        self.Phi = Phi

        self.L = daeDomain("L",  self, um, "Characteristic dimension (size) of crystals")
        self.ni = daeVariable("ni", pbm_number_density_t, self, "Van Leer k-interpolation scheme (k = 1/3)", [self.L])

    def DeclareEquations(self):
        daeModel.DeclareEquations(self)

        G  = self.G 
        hr = daeHRUpwindSchemeEquation(self.ni, self.L, self.Phi, Constant(1e-10 * pbm_number_density_t.Units))
        L  = self.L.Points
        nL = self.L.NumberOfPoints

        # High-resolution cell-centered finite-volume upwind scheme.
        # Available functions:
        #  - dc_dt:   accumulation term (cell-average quantity assumed)
        #  - dc_dx:   convection term (can include the source term integral too)
        #  - d2c_dx2: diffusion term
        for i in range(1, nL):
            eq = self.CreateEquation("ni(%d)" % i, "")
            eq.Residual = hr.dc_dt(i) + G * hr.dc_dx(i, S = None)
                
        # Boundary condition: here G*ni = 0
        eq = self.CreateEquation("ni(0)")
        eq.Residual = self.ni(0) 

class simBatchReactor(daeSimulation):
    def __init__(self, modelName, N, L, G, ni_0, Phi):
        daeSimulation.__init__(self)
        self.m = modelMoC(modelName, G, Phi)
        self.m.Description = __doc__

        self.N    = N
        self.L    = L
        self.ni_0 = ni_0

    def SetUpParametersAndDomains(self):
        self.m.L.CreateStructuredGrid(self.N, min(self.L), max(self.L))
        self.m.L.Points = self.L

    def SetUpVariables(self):
        # Initial conditions
        L = self.m.L.Points
        nL = self.m.L.NumberOfPoints

        # Initial conditions (the first item is not differential)
        self.ni_0[0] = None
        self.m.ni.SetInitialConditions(self.ni_0)

def run_simulation(simPrefix, modelName, N, L, G, ni_0, reportingInterval, timeHorizon, Phi):
    # Create Log, Solver, DataReporter and Simulation object
    log          = daePythonStdOutLog()
    daesolver    = daeIDAS()
    datareporter = daeDelegateDataReporter()
    dr_tcpip     = daeTCPIPDataReporter()
    dr_data      = daeNoOpDataReporter()
    datareporter.AddDataReporter(dr_tcpip)
    datareporter.AddDataReporter(dr_data)
    simulation   = simBatchReactor(modelName, N, L, G, ni_0, Phi)

    # Enable reporting of all variables
    simulation.m.SetReportingOn(True)

    # Set the time horizon and the reporting interval
    simulation.ReportingInterval = reportingInterval
    simulation.TimeHorizon       = timeHorizon

    # Connect data reporter
    simName = simPrefix + strftime(" [%d.%m.%Y %H:%M:%S]", localtime())
    if(dr_tcpip.Connect("", simName) == False):
        sys.exit()

    # Initialize the simulation
    simulation.Initialize(daesolver, datareporter, log)

    # Save the model report and the runtime model report
    #simulation.m.SaveModelReport(simulation.m.Name + ".xml")
    #simulation.m.SaveRuntimeModelReport(simulation.m.Name + "-rt.xml")

    # Solve at time=0 (initialization)
    simulation.SolveInitial()

    # Run
    simulation.Run()
    simulation.Finalize()
    return dr_data.Process

def run(**kwargs):
    # Create an initial CSD array and growth rate
    N    = 200
    L_lb = 0.0
    L_ub = 100.0
    G    = 0.1 * um/s
    ni_0 = numpy.zeros(N+1)
    L    = numpy.linspace(L_lb, L_ub, N+1)
    sigma = 0.778 * (L[1] - L[0])
    for i in range(0, N+1):
        if L[i] > 0 and L[i] < 2:
            ni_0[i] = 0
        elif L[i] > 2 and L[i] <= 10:
            ni_0[i] = 1e10
        elif L[i] > 10 and L[i] <= 18:
            ni_0[i] = 0
        elif L[i] > 18 and L[i] <= 34:
            ni_0[i] = 1e10 * numpy.cos(numpy.pi * (L[i] - 26)/16)**2
        elif L[i] > 34 and L[i] <= 42:
            ni_0[i] = 0
        elif L[i] > 42 and L[i] <= 58:
            ni_0[i] = 1e10 * numpy.sqrt(1 - ((L[i] - 50)**2)/64)
        elif L[i] > 58 and L[i] <= 66:
            ni_0[i] = 0
        elif L[i] > 66 and L[i] <= 77:
            ni_0[i] = 1e10 * numpy.exp(-((L[i] - 70)**2) / (2 * (sigma**2)))
        else:
            ni_0[i] = 0

    reportingInterval = 5
    timeHorizon       = 100
    timeIndex         = 20

    # First find analytical solution
    process_analytical = run_analytical('Analytical', 'Analytical', N, L, G, ni_0, reportingInterval, timeHorizon)
    ni_analytical = process_analytical.dictVariables['Analytical.ni']

    L_report = []
    try:
        for flux_limiter in daeHRUpwindSchemeEquation.supported_flux_limiters:
            flux_limiter_name = flux_limiter.__doc__
            process_dpb = run_simulation(flux_limiter_name, flux_limiter_name, N, L, G, ni_0, reportingInterval, timeHorizon, flux_limiter)
            ni_dpb      = process_dpb.dictVariables['%s.ni' % flux_limiter_name]

            ni_anal = ni_analytical.Values[timeIndex]
            ni_dpb  = ni_dpb.Values[timeIndex]

            L1 = numpy.linalg.norm(ni_dpb-ni_anal, ord = 1)
            L2 = numpy.linalg.norm(ni_dpb-ni_anal, ord = 2)
            print('L1 = %e, L2 = %e' % (L1, L2))

            L_report.append((flux_limiter_name, L1, L2))

    finally:
        # Sort by L2
        L_report.sort(key = lambda t: t[2])
        print('   -------------------------------------')
        print('           Scheme  L1         L2        ')
        print('   -------------------------------------')
        for flux_limiter, L1, L2 in L_report:
            print('  %15s  %.3e  %.3e' % (flux_limiter, L1, L2))
        print('   -------------------------------------')

if __name__ == "__main__":
    run()