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

"""
***********************************************************************************
                           tutorial_opencs_dae_7_cv.py
                DAE Tools: pyOpenCS 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.

This problem and its solution in COMSOL Multiphysics <https://www.comsol.com> software
is described in the COMSOL blog: Verify Simulations with the Method of Manufactured Solutions (2015),
https://www.comsol.com/blogs/verify-simulations-with-the-method-of-manufactured-solutions.

Here, a 1D transient heat conduction problem in a bar of length L is solved using the
Centered Finite Difference method:
   dT/dt - k/(rho*cp) * d2T/dx2 = 0, x in [0,L]
Dirichlet boundary conditions are applied at both bar ends (500 K).
Initial conditions are 500 K.

The manufactured solution is given by function u(x):
   u(x) = 500 + (x/L) * (x/L - 1) * (t/tau)

Note:
 csSimulate cannot be called multiple times from python (multiple MPI init and finalize are not allowed).
Therefore, the simulation are performed using the csSimulator binary
(it has to be in the PATH for this to work).
"""

import os, sys, json, itertools, numpy, pandas
from daetools.solvers.opencs import pyOpenCS
from daetools.solvers.opencs import csModelBuilder_t, csNumber_t, csSimulate
from daetools.solvers.opencs import csGraphPartitioner_t, createGraphPartitioner_2D_Npde
from daetools.examples.tutorial_opencs_aux import compareResults
import matplotlib.pyplot as plt

class HeatConductionBar_1D:
    def __init__(self, Nx_):
        self.Nx = Nx_ + 1
        Nx = self.Nx

        self.L = 100.0 # m

        self.x0 = 0.0
        self.x1 = self.L
        self.dx = (self.x1-self.x0) / (Nx-1)

        self.x_domain = []
        for x in range(self.Nx):
            self.x_domain.append(self.x0 + x * self.dx)

        self.T_start_index = 0*Nx
        self.u_start_index = 1*Nx

        self.Nequations = 2*Nx

    def GetInitialConditions(self):
        # Use numpy array so that setting T_0 and u_0 changes the original values
        uv0 = numpy.zeros(self.Nequations)
        T_0 = uv0[self.T_start_index : self.u_start_index]
        u_0 = uv0[self.u_start_index : self.Nequations]

        for ix in range(self.Nx):
            xp = self.x_domain[ix]

            T_0[ix] = 500.0
            u_0[ix] = 500.0

        return uv0.tolist()

    def GetVariableNames(self):
        return ['T(%d)'%x for x in range(self.Nx)] + ['u(%d)'%x for x in range(self.Nx)]

    def CreateEquations(self, y, dydt, currentTime):
        # y is a list of csNumber_t objects representing model variables
        # dydt is a list of csNumber_t objects representing time derivatives of model variables
        T_values    = y   [self.T_start_index : self.u_start_index]
        u_values    = y   [self.u_start_index : self.Nequations]
        dTdt_values = dydt[self.T_start_index : self.u_start_index]
        
        Nx = self.Nx
        dx = self.dx

        #Ac    =    0.1  # m**2
        rho   = 2700.0  # kg/m**3
        cp    =  900.0  # J/(kg*K)
        kappa =  238.0  # W/(m*K)
        tau   = 3600.0  # seconds
        L     =  self.L # m
        t     = currentTime
        # Thermal diffusivity (m**2/s)
        alpha = kappa/(rho * cp)

        # Manufactured solution
        # Note: x is not an index here, but length in meters
        u       = lambda x: 500.0 + (x/L) * (x/L - 1) * (t/tau)
        du_dt   = lambda x: (x/L) * (x/L - 1) * (1/tau)
        d2u_dx2 = lambda x: (2/(L**2)) * (t/tau)
        Q       = lambda x: du_dt(x) - alpha * d2u_dx2(x)

        # Variables
        T       = lambda x:    T_values[x]
        um      = lambda x:    u_values[x]
        dT_dt   = lambda x: dTdt_values[x]

        def d2T_dx2(x):
            # Interior points
            Ti1 = T(x+1)
            Ti  = T(x)
            Ti2 = T(x-1)
            return (Ti1 - 2*Ti + Ti2) / (dx*dx)

        eq = 0
        equations = [None] * self.Nequations
        # Component T:
        for x in range(Nx):
            if(x == 0):       # BCs: Dirichlet
                equations[eq] = T(x) - 500.0
            elif(x == Nx-1):  # BCs: Dirichlet
                equations[eq] = T(x) - 500.0
            else:             # Interior points
                xp = self.x_domain[x]
                equations[eq] = dT_dt(x) - alpha * d2T_dx2(x) - Q(xp)
            eq += 1
        
        # Manuf. solution u:
        for x in range(Nx):
            xp = self.x_domain[x]
            equations[eq] = um(x) - u(xp)
            eq += 1
        
        return equations

def loadOpenCSData(Nx, resultsDirectory = 'tutorial_opencs_dae_7_cv'):
    csvFilename = os.path.join(resultsDirectory, '%d' % Nx, 'results-0.csv')
    csv_filepath = os.path.join(os.path.dirname(os.path.abspath(__file__)), csvFilename)
    df = pandas.read_csv(csv_filepath, sep=';', header=2, skiprows=None, quotechar='"', skipinitialspace=True, dtype=float)

    Nx = Nx+1

    print('Loading %d data...' % (Nx-1))

    T_start = 1 # The first value is Time (skip it)
    T_end   = T_start + Nx

    u_start = T_end
    u_end   = u_start + Nx

    data = df.values[-1] # take the last time point
    T  =  data[T_start:T_end]
    um =  data[u_start:u_end]

    return T, um

def simulate(Nx):
    # Instantiate the model being simulated.
    model = HeatConductionBar_1D(Nx)
    
    # 1. Initialise the DAE system with the number of variables and other inputs.
    mb = csModelBuilder_t()
    mb.Initialize_DAE_System(model.Nequations, 0, defaultAbsoluteTolerance = 1e-6)
    
    # 2. Specify the OpenCS model.
    mb.ModelEquations = model.CreateEquations(mb.Variables, mb.TimeDerivatives, mb.Time)
    # Set initial conditions.
    mb.VariableValues = model.GetInitialConditions()
    # Set variable names.
    mb.VariableNames  = model.GetVariableNames()
    
    # 3. Generate a model for single CPU simulations.    
    # Set simulation options (here, specified as a dictionary).
    options = mb.SimulationOptions
    options['Simulation']['OutputDirectory']             = '%d' % Nx
    options['Simulation']['TimeHorizon']                 = 20*3600.0
    options['Simulation']['ReportingInterval']           = 3600.0
    options['Solver']['Parameters']['RelativeTolerance'] = 1e-5
    
    options['LinearSolver'] = {
        "Library":   "Trilinos",
        "Name":      "Amesos_Klu",
        "PrintInfo": False,
        "Parameters": {}
    }
    mb.SimulationOptions = options
    
    # Partition the system to create the OpenCS model for a single CPU simulation.
    # In this case (Npe = 1) the graph partitioner is not required.
    Npe = 1
    graphPartitioner = None
    inputFilesDirectory = 'tutorial_opencs_dae_7_cv'
    cs_models = mb.PartitionSystem(Npe, graphPartitioner)
    csModelBuilder_t.ExportModels(cs_models, inputFilesDirectory, mb.SimulationOptions)
    print("Single CPU OpenCS model generated successfully!")

    # 4. Run simulation using the exported model from the specified directory.
    # Cannot call csSimulate multiple times from python (multiple MPI init and finalize are not allowed).
    # Therefore, simulate it using the csSimulator binary (it has to be in the PATH for this to work).
    #csSimulate(inputFilesDirectory)
    os.system('csSimulator %s' % inputFilesDirectory)

    x_coords = numpy.linspace(0, model.L, Nx+1)
    T, um = loadOpenCSData(Nx)
    return  x_coords, T, um

def run(**kwargs):
    Nx1 = 5
    Nx2 = 20
    L = 100.0
    h1 = L / Nx1
    h2 = L / Nx2
    xcoords1, T1, u1 = simulate(Nx1)
    xcoords2, T2, u2 = simulate(Nx2)

    # The normalised global errors
    E1 = numpy.sqrt((1.0/Nx1) * numpy.sum((T1-u1)**2))
    E2 = numpy.sqrt((1.0/Nx2) * numpy.sum((T2-u2)**2))

    # Order of accuracy
    p = numpy.log(E1/E2) / numpy.log(h1/h2)
    C = E1 / h1**p

    print('\n\nOrder of Accuracy:')
    print('||E(h)|| is proportional to: C * (h**p)')
    print('||E(h1)|| = %e, ||E(h2)|| = %e' % (E1, E2))
    print('C = %e' % C)
    print('Order of accuracy (p) = %.2f' % p)

    fontsize = 14
    fontsize_legend = 11
    fig = plt.figure(figsize=(10,4), facecolor='white')
    #title = 'Plots for coarse and fine grids (Order of accuracy = %.2f) (cv_5)' % (p)
    #fig.canvas.setWindowTitle(title)

    ax = plt.subplot(121)
    plt.plot(xcoords1, T1, 'rs', label='T (CFDM)')
    plt.plot(xcoords1, u1, 'b-', label='u (manufactured)')
    plt.xlabel('x, m', fontsize=fontsize)
    plt.ylabel('Temperature, K', fontsize=fontsize)
    plt.legend(fontsize=fontsize_legend)
    plt.xlim((0, 100))

    ax = plt.subplot(122)
    plt.plot(xcoords2, T2, 'rs', label='T (CFDM)')
    plt.plot(xcoords2, u2, 'b-', label='u (manufactured)')
    plt.xlabel('x, m', fontsize=fontsize)
    plt.ylabel('Temperature, K', fontsize=fontsize)
    plt.legend(fontsize=fontsize_legend)
    plt.xlim((0, 100))

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    run()