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

"""
***********************************************************************************
                           tutorial_opencs_dae_8_kernels.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__ = """
In this example the Cahn-Hilliard equation is solved using the finite difference method.
This equation describes the process of phase separation, where two components of a
binary mixture separate and form domains pure in each component.

   dc/dt = Diffusivity * nabla^2(mu)
   mu = c^3 - c - gamma * nabla^2(c)

The mesh is a simple square (0-100)x(0-100).
Input parameters are Diffusivity = 1 and gamma = 1.
For both c an mu insulated boundary conditions are set (no flux on boundaries).
Initial conditions are set to c(0) = 0.0 + c_noise where the noise is specified using
the normal distribution with standard deviation of 0.1.
The system is integrated for 500 seconds and the outputs are taken every 5 seconds.

This version is implemented using kernels.
"""

import os, sys, itertools, numpy, pandas, time
from daetools.solvers.opencs import csModelBuilder_t, csSimulate
from daetools.solvers.opencs import csGroup_t, csKernel_t, csEquation_t
import matplotlib.pyplot as plt
from matplotlib import animation

Diffusivity = 1.0
gamma       = 1.0

class CahnHilliard_Kernels_2D:
    def __init__(self, Nx, Ny, c_flux_bc, mu_flux_bc):
        self.Nx = Nx
        self.Ny = Ny
        self.c_flux_bc  = c_flux_bc
        self.mu_flux_bc = mu_flux_bc
        
        self.x0 = 0.0
        self.x1 = 100.0
        self.y0 = 0.0
        self.y1 = 100.0
        self.dx = (self.x1-self.x0) / (Nx-1)
        self.dy = (self.y1-self.y0) / (Ny-1)

        self.x_domain = []
        self.y_domain = []
        for x in range(self.Nx):
            self.x_domain.append(self.x0 + x * self.dx)
        for y in range(self.Ny):
            self.y_domain.append(self.y0 + y * self.dy)

        self.c_start_index  = 0*Nx*Ny
        self.mu_start_index = 1*Nx*Ny

        self.Nequations = 2*Nx*Ny

        self.group_BCs = csGroup_t ("BoundaryConditions", 1) # the group for boundary conditions (group 1)
        self.kernel_c  = csKernel_t("CahnHilliard_c",     2) # the kernel for c-component (group 2)
        self.kernel_mu = csKernel_t("CahnHilliard_mu",    3) # the kernel for mu-component (group 3)

    def GetInitialConditions(self):
        # Use numpy array so that setting u_0 and v_0 changes the original values
        uv0  = numpy.zeros(self.Nequations)
        c_0  = uv0[self.c_start_index  : self.mu_start_index]
        mu_0 = uv0[self.mu_start_index : self.Nequations]

        numpy.random.seed(124)
        c0     = 0.0
        stddev = 0.1
        def c_with_noise():
            return numpy.random.normal(c0, stddev)

        for ix in range(self.Nx):
            for iy in range(self.Ny):
                index = self.GetIndex(ix,iy)

                c_0[index] = c_with_noise()

        return uv0.tolist()

    def GetVariableNames(self):
        x_y_inds = [(x,y) for x,y in itertools.product(range(self.Nx), range(self.Ny))]
        return ['c(%d,%d)'%(x,y) for x,y in x_y_inds] + ['mu(%d,%d)'%(x,y) for x,y in x_y_inds]

    def CreateEquations(self, y, dydt):
        # 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
        c_values     = y   [self.c_start_index  : self.mu_start_index]
        mu_values    = y   [self.mu_start_index : self.Nequations]
        dcdt_values  = dydt[self.c_start_index  : self.mu_start_index]
        dmudt_values = dydt[self.mu_start_index : self.Nequations]
        
        c_flux_bc  = self.c_flux_bc
        mu_flux_bc = self.mu_flux_bc
        
        dx = self.dx
        dy = self.dy
        Nx = self.Nx
        Ny = self.Ny

        def c(x, y):
            index = self.GetIndex(x, y)
            return c_values[index]
        
        def mu(x, y):
            index = self.GetIndex(x, y)
            return mu_values[index]

        def dc_dt(x, y):
            index = self.GetIndex(x, y)
            return dcdt_values[index]

        def dmu_dt(x, y):
            index = self.GetIndex(x, y)
            return dmudt_values[index]

        # First order partial derivative per x.
        def dc_dx(x, y):
            if(x == 0): # left
                u0 = c(0, y)
                u1 = c(1, y)
                u2 = c(2, y)
                return (-3*u0 + 4*u1 - u2) / (2*dx)
            elif(x == Nx-1): # right
                un  = c(Nx-1,   y)
                un1 = c(Nx-1-1, y)
                un2 = c(Nx-1-2, y)
                return (3*un - 4*un1 + un2) / (2*dx)
            else:
                u1 = c(x+1, y)
                u2 = c(x-1, y)
                return (u1 - u2) / (2*dx)

        def dmu_dx(x, y):
            if(x == 0): # left
                u0 = mu(0, y)
                u1 = mu(1, y)
                u2 = mu(2, y)
                return (-3*u0 + 4*u1 - u2) / (2*dx)
            elif(x == Nx-1): # right
                un  = mu(Nx-1,   y)
                un1 = mu(Nx-1-1, y)
                un2 = mu(Nx-1-2, y)
                return (3*un - 4*un1 + un2) / (2*dx)
            else:
                u1 = mu(x+1, y)
                u2 = mu(x-1, y)
                return (u1 - u2) / (2*dx)

        # First order partial derivative per y.
        def dc_dy(x, y):
            if(y == 0): # bottom
                u0 = c(x, 0)
                u1 = c(x, 1)
                u2 = c(x, 2)
                return (-3*u0 + 4*u1 - u2) / (2*dy)
            elif(y == Ny-1): # top
                un  = c(x, Ny-1  )
                un1 = c(x, Ny-1-1)
                un2 = c(x, Ny-1-2)
                return (3*un - 4*un1 + un2) / (2*dy)
            else:
                ui1 = c(x, y+1)
                ui2 = c(x, y-1)
                return (ui1 - ui2) / (2*dy)

        def dmu_dy(x, y):
            if(y == 0): # bottom
                u0 = mu(x, 0)
                u1 = mu(x, 1)
                u2 = mu(x, 2)
                return (-3*u0 + 4*u1 - u2) / (2*dy)
            elif(y == Ny-1): # top
                un  = mu(x, Ny-1  )
                un1 = mu(x, Ny-1-1)
                un2 = mu(x, Ny-1-2)
                return (3*un - 4*un1 + un2) / (2*dy)
            else:
                ui1 = mu(x, y+1)
                ui2 = mu(x, y-1)
                return (ui1 - ui2) / (2*dy)

        # Second order partial derivative per x.
        def d2c_dx2(x, y):
            if(x == 0 or x == Nx-1):
                raise RuntimeError("d2c_dx2 called at the boundary")

            ui1 = c(x+1, y)
            ui  = c(x,   y)
            ui2 = c(x-1, y)
            return (ui1 - 2*ui + ui2) / (dx*dx)

        def d2mu_dx2(x, y):
            if(x == 0 or x == Nx-1):
                raise RuntimeError("d2mu_dx2 called at the boundary")

            vi1 = mu(x+1, y)
            vi  = mu(x,   y)
            vi2 = mu(x-1, y)
            return (vi1 - 2*vi + vi2) / (dx*dx)

        # Second order partial derivative per y.
        def d2c_dy2(x, y):
            if(y == 0 or y == Ny-1):
                raise RuntimeError("d2c_dy2 called at the boundary")

            ui1 = c(x, y+1)
            ui  = c(x,   y)
            ui2 = c(x, y-1)
            return (ui1 - 2*ui + ui2) / (dy*dy)

        def d2mu_dy2(x, y):
            if(y == 0 or y == Ny-1):
                raise RuntimeError("d2mu_dy2 called at the boundary")

            vi1 = mu(x, y+1)
            vi  = mu(x,   y)
            vi2 = mu(x, y-1)
            return (vi1 - 2*vi + vi2) / (dy*dy)

        equations = []
        kernels   = [self.kernel_c, self.kernel_mu]
        # Component c:
        for x in range(Nx):
            for y in range(Ny):
                if(x == 0):       # Left BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ c(x,y) ] = dc_dx(x,y) - c_flux_bc
                    equations.append(equation)
                elif(x == Nx-1):  # Right BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ c(x,y) ] = dc_dx(x,y) - c_flux_bc
                    equations.append(equation)
                elif(y == 0):     # Bottom BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ c(x,y) ] = dc_dy(x,y) - c_flux_bc
                    equations.append(equation)
                elif(y == Ny-1):  # Top BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ c(x,y) ] = dc_dy(x,y) - c_flux_bc
                    equations.append(equation)
                else:
                    # Interior points
                    equation_c = csEquation_t(self.kernel_c)
                    equation_c[ c(x,y) ] = dc_dt(x,y) - Diffusivity * (d2mu_dx2(x,y) + d2mu_dy2(x,y))
                    equation_c.SetGridPoint(x, y);
                    self.kernel_c.AddEquation(equation_c);
        
        # Component mu:
        for x in range(Nx):
            for y in range(Ny):
                if(x == 0):       # Left BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ mu(x,y) ] = dmu_dx(x,y) - mu_flux_bc
                    equations.append(equation)
                elif(x == Nx-1):  # Right BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ mu(x,y) ] = dmu_dx(x,y) - mu_flux_bc
                    equations.append(equation)
                elif(y == 0):     # Bottom BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ mu(x,y) ] = dmu_dy(x,y) - mu_flux_bc
                    equations.append(equation)
                elif(y == Ny-1):  # Top BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ mu(x,y) ] = dmu_dy(x,y) - mu_flux_bc
                    equations.append(equation)
                else:
                    # Interior points
                    equation_mu = csEquation_t(self.kernel_mu)
                    equation_mu[ mu(x,y) ] = mu(x,y) \
                                             + gamma * (d2c_dx2(x,y) + d2c_dy2(x,y)) \
                                             - (c(x,y)*c(x,y)*c(x,y) - c(x,y))
                    equation_mu.SetGridPoint(x, y);
                    self.kernel_mu.AddEquation(equation_mu);

        return (equations, kernels)
    
    def GetIndex(self, x, y):
        if x < 0 or x >= self.Nx:
            raise RuntimeError("Invalid x index")
        if y < 0 or y >= self.Ny:
            raise RuntimeError("Invalid y index")
        return self.Ny*x + y
    
def run(**kwargs):
    inputFilesDirectory = kwargs.get('inputFilesDirectory', os.path.splitext(os.path.basename(__file__))[0])
    Nx         = kwargs.get('Nx',         100)
    Ny         = kwargs.get('Ny',         100)
    c_flux_bc  = kwargs.get('c_flux_bc',  0.0)
    mu_flux_bc = kwargs.get('mu_flux_bc', 0.0)
    
    # Instantiate the model being simulated.
    model = CahnHilliard_Kernels_2D(Nx, Ny, c_flux_bc, mu_flux_bc)
    
    # 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-5)
    
    # 2. Specify the OpenCS model.
    # Create and set model equations using the provided time/variable/timeDerivative/dof objects.
    # The DAE system is defined as:
    #     F(x',x,y,t) = 0
    # where x' are derivatives of state variables, x are state variables,
    # y are fixed variables (degrees of freedom) and t is the current simulation time.
    mb.ModelEquations = model.CreateEquations(mb.Variables, mb.TimeDerivatives)    
    # 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 (specified as a string in JSON format).
    options = mb.SimulationOptions
    options['Simulation']['OutputDirectory']             = 'results'
    options['Simulation']['TimeHorizon']                 = 500.0
    options['Simulation']['ReportingInterval']           =   5.0
    options['Solver']['Parameters']['RelativeTolerance'] =  1e-5
    
    # (1) Sequential evaluator (evaluates groups one by one).
    options['Model']['Evaluators']['Device_0']['Library'] = 'Sequential'
    options['Model']['Evaluators']['Device_0']['Name']    = 'seq'
    options['Model']['Evaluators']['Device_0']['Groups']  = ['BoundaryConditions']
    # (2a) OpenMP evaluator.
    options['Model']['Evaluators']['Device_1']               = {}
    options['Model']['Evaluators']['Device_1']['Parameters'] = {}
    options['Model']['Evaluators']['Device_1']['Library'] = 'Kernels_OpenMP'
    options['Model']['Evaluators']['Device_1']['Name']    = 'kopenmp'
    options['Model']['Evaluators']['Device_1']['Groups']  = []
    options['Model']['Evaluators']['Device_1']['Kernels'] = ['CahnHilliard_c', 'CahnHilliard_mu']
    options['Model']['Evaluators']['Device_1']['Parameters']['numThreads'] = 0

    # ILU options for Ncpu = 1: k = 3, rho = 1.0, alpha = 1e-1, w = 0.5
    options['LinearSolver']['Preconditioner']['Parameters']['fact: level-of-fill']      =    3
    options['LinearSolver']['Preconditioner']['Parameters']['fact: relax value']        =  0.0
    options['LinearSolver']['Preconditioner']['Parameters']['fact: absolute threshold'] = 1e-1
    options['LinearSolver']['Preconditioner']['Parameters']['fact: relative threshold'] =  1.0
    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
    cs_models = mb.PartitionSystem(Npe, graphPartitioner)
    csModelBuilder_t.ExportModels(cs_models, inputFilesDirectory, mb.SimulationOptions)
    print("Single CPU OpenCS model generated successfully!")

    # 5. Run simulation using the exported model from the specified directory.
    csSimulate(inputFilesDirectory)
    
    try:
        csvFilename = os.path.join(inputFilesDirectory, 'results', '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)

        # The data layout in memory is as (for each time point):
        #   [Time][BoundaryConditions][CahnHilliard_c][CahnHilliard_mu]
        #   [  1 ][  2*(2*Nx+2*Ny-4) ][ (Nx-2)*(Ny-2)][ (Nx-2)*(Ny-2) ]

        # Offset the data by time and the whole BoundaryConditions group.
        c_start   = 1 + 2*(2*Nx + 2*Ny - 4)
        c_end     = c_start + (Nx-2)*(Ny-2)

        mu_start   = c_end
        mu_end     = mu_start + (Nx-2)*(Ny-2)

        # a) For simplicity, plot it without the boundaries.
        def getData_c_no_bc(data):
            c = data[c_start : c_end]
            c = numpy.array(c).reshape(Nx-2,Ny-2)
            return c

        # b) Add both boundary and the interior points.
        def getData_c(data):
            cnt = 1 # The first point is always the Time.
            c = numpy.zeros((Nx,Ny), dtype=float)
            # The equations are sorted by groups so that the boundary points are located
            # at the beggining of the data array.
            # First, iterate over boundary points and set their values
            for x in range(Nx):
                for y in range(Ny):
                    if (x == 0) or (x == Nx-1) or (y == 0) or (y == Ny-1):
                        c[x,y] = data[cnt]
                        cnt += 1
                    else:
                        continue
            # Second, set the interior points (the group CahnHilliard_c is after the BoundaryConditions)
            c[1:Nx-1,1:Ny-1] = numpy.array(data[c_start : c_end]).reshape(Nx-2,Ny-2)
            return c

        Nframes = len(df.values)

        # Plot setup and the first frame
        data = df.values[0]
        t = data[0]
        c = getData_c(data)
        fig = plt.figure()
        im = plt.imshow(c, cmap = 'viridis')
        plt.colorbar()
        plt.title('time %.2f s' % t)

        def animate(i):
            data = df.values[i]
            t = data[0]
            c = getData_c(data)
            im.set_data(c)
            plt.title('time %.2f s' % t)
            return im

        # Save the animation object so it does not go out if scope while animation is running.
        plt.__anim__ = animation.FuncAnimation(fig, animate, frames=Nframes, interval = 500, repeat = False)
        plt.show()
    except Exception as e:
        print(str(e))

if __name__ == "__main__":
    if len(sys.argv) == 1:
        Nx = 100
        Ny = 100
    elif len(sys.argv) == 3:
        Nx = int(sys.argv[1])
        Ny = int(sys.argv[2])
    else:
        print('Usage: python tutorial_opencs_dae_8_kernels.py Nx Ny')
        sys.exit()
        
    c_flux_bc  = 0.0
    mu_flux_bc = 0.0
    inputFilesDirectory = 'tutorial_opencs_dae_8_kernels'
    run(Nx = Nx, 
        Ny = Ny, 
        c_flux_bc = c_flux_bc,
        mu_flux_bc = mu_flux_bc,
        inputFilesDirectory = inputFilesDirectory)