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

"""
***********************************************************************************
                           tutorial_opencs_dae_3_groups.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__ = """
Reimplementation of tutorial_opencs_dae_3 using groups.
"""

import os, sys, json, itertools, numpy
from daetools.solvers.opencs import csModelBuilder_t, csNumber_t, csSimulate
from daetools.solvers.opencs import csGroup_t, csKernel_t, csEquation_t
from daetools.solvers.opencs import csGraphPartitioner_t, createGraphPartitioner_2D_Npde
from daetools.examples.tutorial_opencs_aux import compareResults

eps1 = 0.002
eps2 = 0.002
A    = 1.000
B    = 3.400

class BrusselatorGroups_2D:
    def __init__(self, Nx, Ny, u_flux_bc, v_flux_bc):
        self.Nx = Nx
        self.Ny = Ny
        self.u_flux_bc = u_flux_bc
        self.v_flux_bc = v_flux_bc
        
        self.x0 = 0.0
        self.x1 = 1.0
        self.y0 = 0.0
        self.y1 = 1.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.u_start_index = 0*Nx*Ny
        self.v_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.group_u   = csGroup_t("Brusselator_u",      2) # the group for u-component (group 2)
        self.group_v   = csGroup_t("Brusselator_v",      3) # the group for v-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)
        u_0 = uv0[self.u_start_index : self.v_start_index]
        v_0 = uv0[self.v_start_index : self.Nequations]

        dx = self.dx
        dy = self.dy
        x0 = self.x0
        x1 = self.x1
        y0 = self.y0
        y1 = self.y1
        
        pi = 3.1415926535898
        Lx = x1 - x0
        Ly = y1 - y0
        
        for ix in range(self.Nx):
            for iy in range(self.Ny):
                index = self.GetIndex(ix,iy)
                x = self.x_domain[ix]
                y = self.y_domain[iy]

                u_0[index] = 1.0 - 0.5 * numpy.cos(pi * y / Ly)
                v_0[index] = 3.5 - 2.5 * numpy.cos(pi * x / Lx)

        return uv0.tolist()

    def GetVariableNames(self):
        x_y_inds = [(x,y) for x,y in itertools.product(range(self.Nx), range(self.Ny))]
        return ['u(%d,%d)'%(x,y) for x,y in x_y_inds] + ['v(%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
        u_values    = y   [self.u_start_index : self.v_start_index]
        v_values    = y   [self.v_start_index : self.Nequations]
        dudt_values = dydt[self.u_start_index : self.v_start_index]
        dvdt_values = dydt[self.v_start_index : self.Nequations]
        
        u_flux_bc = self.u_flux_bc
        v_flux_bc = self.v_flux_bc
        
        dx = self.dx
        dy = self.dy
        Nx = self.Nx
        Ny = self.Ny
        x_domain = self.x_domain
        y_domain = self.y_domain

        def u(x, y):
            index = self.GetIndex(x, y)
            return u_values[index]
        
        def v(x, y):
            index = self.GetIndex(x, y)
            return v_values[index]

        def du_dt(x, y):
            index = self.GetIndex(x, y)
            return dudt_values[index]

        def dv_dt(x, y):
            index = self.GetIndex(x, y)
            return dvdt_values[index]

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

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

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

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

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

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

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

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

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

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

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

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

        equations = []
        kernels   = []
        # Component u:
        for x in range(Nx):
            for y in range(Ny):
                if(x == 0):       # Left BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ u(x,y) ] = du_dx(x,y) - u_flux_bc
                    equations.append(equation)
                elif(x == Nx-1):  # Right BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ u(x,y) ] = du_dx(x,y) - u_flux_bc
                    equations.append(equation)
                elif(y == 0):     # Bottom BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ u(x,y) ] = du_dy(x,y) - u_flux_bc
                    equations.append(equation)
                elif(y == Ny-1):  # Top BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ u(x,y) ] = du_dy(x,y) - u_flux_bc
                    equations.append(equation)
                else:
                    # Interior points
                    equation_u = csEquation_t(self.group_u)
                    equation_u[ u(x,y) ] = du_dt(x,y) \
                                           - eps1 * (d2u_dx2(x,y) + d2u_dy2(x,y)) \
                                           - (u(x,y)*u(x,y)*v(x,y) - (B+1)*u(x,y) + A)
                    equations.append(equation_u)
        
        # Component v:
        for x in range(Nx):
            for y in range(Ny):
                if(x == 0):       # Left BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ v(x,y) ] = dv_dx(x,y) - v_flux_bc
                    equations.append(equation)
                elif(x == Nx-1):  # Right BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ v(x,y) ] = dv_dx(x,y) - v_flux_bc
                    equations.append(equation)
                elif(y == 0):     # Bottom BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ v(x,y) ] = dv_dy(x,y) - v_flux_bc
                    equations.append(equation)
                elif(y == Ny-1):  # Top BC: Neumann BCs
                    equation = csEquation_t(self.group_BCs)
                    equation[ v(x,y) ] = dv_dy(x,y) - v_flux_bc
                    equations.append(equation)
                else:
                    # Interior points
                    equation_v = csEquation_t(self.group_v)
                    equation_v[ v(x,y) ] = dv_dt(x,y) \
                                           - eps2 * (d2v_dx2(x,y) + d2v_dy2(x,y)) \
                                           + u(x,y)*u(x,y)*v(x,y) - B*u(x,y)
                    equations.append(equation_v)

        #for i, eq in enumerate(equations):
        #    print('%d - %d' % (i, eq.targetVariableIndex))
        #print(len(equations))
        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
    
    
class myGraphPartitioner(csGraphPartitioner_t):
    # Divides vertices into Npe partitions without an analysis 
    def __init__(self):
        csGraphPartitioner_t.__init__(self)

    def GetName(self):
        return 'myGraphPartitioner'
    
    def Partition(self, Npe, Nvertices, Nconstraints, rowIndices, colIndices, vertexWeights):
        if Npe > Nvertices:
            raise RuntimeError('Npe larger than Nvertices')
        
        Neq_pe     = int(Nvertices / Npe)
        equations  = list(range(Nvertices))
        partitions = []
        for i in range(Npe):
            start = i     * Neq_pe
            end   = (i+1) * Neq_pe
            if i == Npe-1:
                end = Nvertices
            partitions.append(equations[start:end])
        print('myGraphPartitioner partitions: %s' % partitions)
        return partitions

def run(**kwargs):
    inputFilesDirectory = kwargs.get('inputFilesDirectory', os.path.splitext(os.path.basename(__file__))[0])
    Nx        = kwargs.get('Nx',        82)
    Ny        = kwargs.get('Ny',        82)
    u_flux_bc = kwargs.get('u_flux_bc', 0.0)
    v_flux_bc = kwargs.get('v_flux_bc', 0.0)
    
    # Instantiate the model being simulated.
    model = BrusselatorGroups_2D(Nx, Ny, u_flux_bc, v_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']                 = 20.0
    options['Simulation']['ReportingInterval']           =  0.5
    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'] = 'OpenMP'
    options['Model']['Evaluators']['Device_1']['Name']    = 'openmp'
    options['Model']['Evaluators']['Device_1']['Groups']  = ['Brusselator_u', 'Brusselator_v']
    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.5
    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!")
    
    # 4. Run simulation using the exported model from the specified directory.
    csSimulate(inputFilesDirectory)
    
    # Compare OpenCS and the original model results.
    compareResults(inputFilesDirectory, ['u(0,0)', 'u(81,81)'])
    
if __name__ == "__main__":
    if len(sys.argv) == 1:
        Nx = 82
        Ny = 82
    elif len(sys.argv) == 3:
        Nx = int(sys.argv[1])
        Ny = int(sys.argv[2])
    else:
        print('Usage: python tutorial_opencs_dae_3_groups.py Nx Ny')
        sys.exit()
        
    u_flux_bc = 0.0
    v_flux_bc = 0.0
    inputFilesDirectory = 'tutorial_opencs_dae_3_groups'
    run(Nx = Nx, 
        Ny = Ny, 
        u_flux_bc = u_flux_bc, 
        v_flux_bc = v_flux_bc, 
        inputFilesDirectory = inputFilesDirectory)