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

"""
***********************************************************************************
                           tutorial_cv_2.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__ = """
Code verification using the Method of Manufactured Solutions.

References: 

1. G. Tryggvason. Method of Manufactured Solutions, Lecture 33: Predictivity-I, 2011.
   `PDF <http://www3.nd.edu/~gtryggva/CFD-Course/2011-Lecture-33.pdf>`_
2. K. Salari and P. Knupp. Code Verification by the Method of Manufactured Solutions. 
   SAND2000 – 1444 (2000).
   `doi:10.2172/759450 <https://doi.org/10.2172/759450>`_
3. P.J. Roache. Fundamentals of Verification and Validation. Hermosa, 2009.
   `ISBN-10:0913478121 <http://www.isbnsearch.org/isbn/0913478121>`_

The procedure for the *transient convection-diffusion* (Burger's) equation:
    
.. code-block:: none

   L(f) = df/dt + f*df/dx - D*d2f/dx2 = 0

is the following:
    
1. Pick a function (q, the manufactured solution): 
    
   .. code-block:: none
    
      q = A + sin(x + Ct)

2. Compute the new source term (g) for the original problem:
    
   .. code-block:: none
      
      g = dq/dt + q*dq/dx - D*d2q/dx2

3. Solve the original problem with the new source term:
    
   .. code-block:: none

      df/dt + f*df/dx - D*d2f/dx2 = g

Since L(f) = g and g = L(q), consequently we have: f = q.
Therefore, the computed numerical solution (f) should be equal to the manufactured one (q).

The terms in the source g term are:

.. code-block:: none

   dq/dt   = C * cos(x + C*t)
   dq/dx   = cos(x + C*t)
   d2q/dx2 = -sin(x + C*t)

The equations are solved for Dirichlet boundary conditions:
    
.. code-block:: none

   f(x=0)   = q(x=0)   = A + sin(0 + Ct)
   f(x=2pi) = q(x=2pi) = A + sin(2pi + Ct)

Numerical vs. manufactured solution plot (no. elements = 60, t = 1.0s):

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

The normalised global errors and the order of accuracy plots (no. elements = [60, 90, 120, 150], t = 1.0s):

.. image:: _static/tutorial_cv_2-results2.png
   :width: 800px
"""

import sys, numpy
from time import localtime, strftime
from daetools.pyDAE import *
import matplotlib.pyplot as plt

no_t = daeVariableType("no_t", dimless, -1.0e+20, 1.0e+20, 0.0, 1e-6)

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

        self.x  = daeDomain("x", self, m, "X axis domain")

        self.A = 1.0
        self.C = 1.0
        self.D = 0.05

        self.f = daeVariable("f", no_t, self, "", [self.x])
        self.q = daeVariable("q", no_t, self, "", [self.x])

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

        # Create some auxiliary functions to make equations more readable 
        A       = self.A
        C       = self.C
        D       = self.D
        t       = Time()
        f       = lambda x:    self.f(x)
        df_dt   = lambda x: dt(self.f(x))
        df_dx   = lambda x:  d(self.f(x), self.x, eCFDM)
        d2f_dx2 = lambda x: d2(self.f(x), self.x, eCFDM)
        q       = lambda x: A + numpy.sin(x() + C*t)
        dq_dt   = lambda x: C * numpy.cos(x() + C*t)
        dq_dx   = lambda x: numpy.cos(x() + C*t)
        d2q_dx2 = lambda x: -numpy.sin(x() + C*t)
        g       = lambda x: dq_dt(x) + q(x) * dq_dx(x) - D * d2q_dx2(x)

        # Numerical solution
        eq = self.CreateEquation("f", "Numerical solution")
        x = eq.DistributeOnDomain(self.x, eOpenOpen)
        eq.Residual = df_dt(x) + f(x) * df_dx(x) - D * d2f_dx2(x) - g(x)
        eq.CheckUnitsConsistency = False

        eq = self.CreateEquation("f(0)", "Numerical solution")
        x = eq.DistributeOnDomain(self.x, eLowerBound)
        eq.Residual = f(x) - q(x)
        eq.CheckUnitsConsistency = False

        eq = self.CreateEquation("f(2pi)", "Numerical solution")
        x = eq.DistributeOnDomain(self.x, eUpperBound)
        eq.Residual = f(x) - q(x)
        eq.CheckUnitsConsistency = False

        # Manufactured solution
        eq = self.CreateEquation("q", "Manufactured solution")
        x = eq.DistributeOnDomain(self.x, eClosedClosed)
        eq.Residual = self.q(x) - q(x)
        eq.CheckUnitsConsistency = False

class simTutorial(daeSimulation):
    def __init__(self, Nx):
        daeSimulation.__init__(self)
        self.m = modTutorial("tutorial_cv_2(%d)" % Nx)
        self.m.Description = __doc__
        
        self.Nx = Nx

    def SetUpParametersAndDomains(self):
        self.m.x.CreateStructuredGrid(self.Nx, 0, 2*numpy.pi)
        
    def SetUpVariables(self):
        Nx = self.m.x.NumberOfPoints
        xp = self.m.x.Points
        for x in range(1, Nx-1):
            self.m.f.SetInitialCondition(x, self.m.A + numpy.sin(xp[x]))
                
# Setup everything manually and run in a console
def simulate(Nx):
    # Create Log, Solver, DataReporter and Simulation object
    log          = daePythonStdOutLog()
    daesolver    = daeIDAS()
    datareporter = daeDelegateDataReporter()
    simulation   = simTutorial(Nx)

    # Do no print progress
    log.PrintProgress = False

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

    # Enable reporting of time derivatives for all reported variables
    simulation.ReportTimeDerivatives = True

    # Set the time horizon and the reporting interval
    simulation.ReportingInterval = 0.05
    simulation.TimeHorizon = 1

    # Connect data reporter
    simName = simulation.m.Name + strftime(" [%d.%m.%Y %H:%M:%S]", localtime())

    # 1. TCP/IP
    tcpipDataReporter = daeTCPIPDataReporter()
    datareporter.AddDataReporter(tcpipDataReporter)
    if not tcpipDataReporter.Connect("", simName):
        sys.exit()

    # 2. Data
    dr = daeNoOpDataReporter()
    datareporter.AddDataReporter(dr)

    # 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()
    
    ###########################################
    #  Data                                   #
    ###########################################
    results = dr.Process.dictVariables
    fvar = results[simulation.m.Name + '.f']
    qvar = results[simulation.m.Name + '.q']
    times = fvar.TimeValues
    q = qvar.Values[-1, :] # 2D array [t,x]
    f = fvar.Values[-1, :] # 2D array [t,x]
    #print(times,f,q)
    
    return times,f,q

def run(**kwargs):
    Nxs = numpy.array([60, 90, 120, 150])
    n = len(Nxs)
    L = 2*numpy.pi
    hs = L / Nxs
    E = numpy.zeros(n)
    C = numpy.zeros(n)
    p = numpy.zeros(n)
    E2 = numpy.zeros(n)
    
    # The normalised global errors
    for i,Nx in enumerate(Nxs):
        times, numerical_sol, manufactured_sol = simulate(int(Nx))
        E[i] = numpy.sqrt((1.0/Nx) * numpy.sum((numerical_sol-manufactured_sol)**2))

    # Order of accuracy
    for i,Nx in enumerate(Nxs):
        p[i] = numpy.log(E[i]/E[i-1]) / numpy.log(hs[i]/hs[i-1])
        C[i] = E[i] / hs[i]**p[i]
        
    C2 = 0.09 # constant for the second order slope line (to get close to the actual line)
    E2 = C2 * hs**2 # E for the second order slope
    
    fontsize = 14
    fontsize_legend = 11
    fig = plt.figure(figsize=(10,4), facecolor='white')
    fig.canvas.set_window_title('The Normalised global errors and the Orders of accuracy (Nelems = %s)' % Nxs.tolist())
    
    ax = plt.subplot(121)
    plt.figure(1, facecolor='white')
    plt.loglog(hs, E,  'ro', label='E(h)')
    plt.loglog(hs, E2, 'b-', label='2nd order slope')
    plt.xlabel('h', fontsize=fontsize)
    plt.ylabel('||E||', fontsize=fontsize)
    plt.legend(fontsize=fontsize_legend)
    plt.xlim((0.04, 0.11))
        
    ax = plt.subplot(122)
    plt.figure(1, facecolor='white')
    plt.semilogx(hs[1:], p[1:],  'rs-', label='Order of Accuracy (p)')
    plt.xlabel('h', fontsize=fontsize)
    plt.ylabel('p', fontsize=fontsize)
    plt.legend(fontsize=fontsize_legend)
    plt.xlim((0.04, 0.075))
    plt.ylim((2.0, 2.1))
    
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    run()