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

"""
***********************************************************************************
                           tutorial_opencs_ode_3.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_ode_3 using groups.
"""

import os, sys, json, itertools, numpy
from daetools.solvers.opencs import pyOpenCS
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

V   =  1.00E-03
Kh  =  4.00E-06
Kv0 =  1.00E-08
q1  =  1.63E-16
q2  =  4.66E-16
C3  =  3.70E+16
a3  = 22.62
a4  =  7.601

class DiurnalKineticsGroups_2D:
    def __init__(self, Nx, Ny):
        self.Nx = Nx
        self.Ny = Ny
        
        self.x0 =  0.0
        self.x1 = 20.0
        self.y0 = 30.0
        self.y1 = 50.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.C1_start_index = 0*Nx*Ny
        self.C2_start_index = 1*Nx*Ny

        self.Nequations = 2*Nx*Ny

        self.group_C1 = csGroup_t("DiurnalKinetics_C1", 1) # the group for C1-component (group 2)
        self.group_C2 = csGroup_t("DiurnalKinetics_C2", 2) # the group for C2-component (group 3)

    def GetInitialConditions(self):
        # Use numpy array so that setting C1_0 and C2_0 changes the original values
        C0 = numpy.zeros(self.Nequations)
        C1_0 = C0[self.C1_start_index : self.C2_start_index]
        C2_0 = C0[self.C2_start_index : self.Nequations]

        dx = self.dx
        dy = self.dy
        x0 = self.x0
        x1 = self.x1
        y0 = self.y0
        y1 = self.y1
        
        def alfa(x):
            xmid = (x1+x0)/2.0
            cx = 0.1 * (x - xmid)
            cx2 = cx*cx
            return 1 - cx2 + 0.5*cx2*cx2
        
        def beta(y):
            ymid = (y1+y0)/2.0
            cy = 0.1 * (y - ymid)
            cy2 = cy*cy
            return 1 - cy2 + 0.5*cy2*cy2
        
        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]

                C1_0[index] = 1E6  * alfa(x) * beta(y)
                C2_0[index] = 1E12 * alfa(x) * beta(y)
        
        return C0.tolist()

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

    def CreateEquations(self, y, time):
        # y is a list of csNumber_t objects representing model variables
        C1_values = y[self.C1_start_index : self.C2_start_index]
        C2_values = y[self.C2_start_index : self.Nequations]
        
        dx = self.dx
        dy = self.dy
        Nx = self.Nx
        Ny = self.Ny
        x_domain = self.x_domain
        y_domain = self.y_domain

        # Math. functions
        sin    = pyOpenCS.sin
        acos   = pyOpenCS.acos
        exp    = pyOpenCS.exp
        cs_max = pyOpenCS.max

        nonzero = csNumber_t(1E-30)
        
        def Kv(y):
            return Kv0 * numpy.exp(y_domain[y] / 5.0)

        def q3(time):
            w = numpy.arccos(-1.0) / 43200.0
            sinwt = cs_max(nonzero, sin(w*time))
            return exp(-a3 / sinwt)

        def q4(time):
            w = numpy.arccos(-1.0) / 43200.0
            sinwt = cs_max(nonzero, sin(w*time))
            return exp(-a4 / sinwt)

        def R1(c1, c2, time):
            return -q1*c1*C3 - q2*c1*c2 + 2*q3(time)*C3 + q4(time)*c2

        def R2(c1, c2, time):
            return q1*c1*C3 - q2*c1*c2 - q4(time)*c2

        def C1(x, y):
            index = self.GetIndex(x, y)
            return C1_values[index]
        
        def C2(x, y):
            index = self.GetIndex(x, y)
            return C2_values[index]

        # First order partial derivative per x.
        def dC1_dx(x, y):
            # If called for x == 0 or x == Nx-1 return 0.0 (zero flux through boundaries).
            if(x == 0 or x == Nx-1):
                return csNumber_t(0.0)

            ci1 = C1(x+1, y)
            ci2 = C1(x-1, y)
            return (ci1 - ci2) / (2*dx)

        def dC2_dx(x, y):
            # If called for x == 0 or x == Nx-1 return 0.0 (zero flux through boundaries).
            if(x == 0 or x == Nx-1):
                return csNumber_t(0.0)

            ci1 = C2(x+1, y)
            ci2 = C2(x-1, y)
            return (ci1 - ci2) / (2*dx)

        # First order partial derivative per y.
        def dC1_dy(x, y):
            # If called for y == 0 or y == Ny-1 return 0.0 (zero flux through boundaries).
            if(y == 0 or y == Ny-1):
                return csNumber_t(0.0)

            ci1 = C1(x, y+1)
            ci2 = C1(x, y-1)
            return (ci1 - ci2) / (2*dy)

        def dC2_dy(x, y):
            # If called for y == 0 or y == Ny-1  return 0.0 (zero flux through boundaries).
            if(y == 0 or y == Ny-1):
                return csNumber_t(0.0)

            ci1 = C2(x, y+1)
            ci2 = C2(x, y-1)
            return (ci1 - ci2) / (2*dy)

        # Second order partial derivative per x.
        def d2C1_dx2(x, y):
            # If called for x == 0 or x == Nx-1 return 0.0 (no diffusion through boundaries).
            if(x == 0 or x == Nx-1):
                return csNumber_t(0.0)

            ci1 = C1(x+1, y)
            ci  = C1(x,   y)
            ci2 = C1(x-1, y)
            return (ci1 - 2*ci + ci2) / (dx*dx)

        def d2C2_dx2(x, y):
            # If called for x == 0 or x == Nx-1 return 0.0 (no diffusion through boundaries).
            if(x == 0 or x == Nx-1):
                return csNumber_t(0.0)

            ci1 = C2(x+1, y)
            ci  = C2(x,   y)
            ci2 = C2(x-1, y)
            return (ci1 - 2*ci + ci2) / (dx*dx)

        # Second order partial derivative per y.
        def d2C1_dy2(x, y):
            # If called for y == 0 or y == Ny-1 return 0.0 (no diffusion through boundaries).
            if(y == 0 or y == Ny-1):
                return csNumber_t(0.0)

            ci1 = C1(x, y+1)
            ci  = C1(x,   y)
            ci2 = C1(x, y-1)
            return (ci1 - 2*ci + ci2) / (dy*dy)

        def d2C2_dy2(x, y):
            # If called for y == 0 or y == Ny-1 return 0.0 (no diffusion through boundaries).
            if(y == 0 or y == Ny-1):
                return csNumber_t(0.0)

            ci1 = C2(x, y+1)
            ci  = C2(x,   y)
            ci2 = C2(x, y-1)
            return (ci1 - 2*ci + ci2) / (dy*dy)

        equations = []
        kernels   = []
        # Component 2 (C1):
        for x in range(Nx):
            for y in range(Ny):
                equation = csEquation_t(self.group_C1)
                equation[ C1(x,y) ] = V  * dC1_dx(x,y) +  \
                                      Kh * d2C1_dx2(x,y) + \
                                      Kv(y) * (0.2 * dC1_dy(x,y) + d2C1_dy2(x,y)) + \
                                      R1(C1(x,y), C2(x,y), time)
                equations.append(equation)

        # Component 2 (C2):
        for x in range(Nx):
            for y in range(Ny):
                equation = csEquation_t(self.group_C2)
                equation[ C2(x,y) ] = V  * dC2_dx(x,y) + \
                                      Kh * d2C2_dx2(x,y) + \
                                      Kv(y) * (0.2 * dC2_dy(x,y) + d2C2_dy2(x,y)) + \
                                      R2(C1(x,y), C2(x,y), time)
                equations.append(equation)

        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', 80)
    Ny = kwargs.get('Ny', 80)
    
    # Instantiate the model being simulated.
    model = DiurnalKineticsGroups_2D(Nx, Ny)
    
    # 1. Initialise the ODE system with the number of variables and other inputs.
    mb = csModelBuilder_t()
    mb.Initialize_ODE_System(model.Nequations, 0, defaultAbsoluteTolerance = 1e-5)
    
    # Create and set model equations using the provided time/variable/dof objects.
    # The ODE system is defined as:
    #     x' = f(x,y,t)
    # 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.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 (specified as a string in JSON format).
    options = mb.SimulationOptions
    options['Simulation']['OutputDirectory']             = 'results'
    options['Simulation']['TimeHorizon']                 = 86400.0
    options['Simulation']['ReportingInterval']           =   360.0
    options['Solver']['Parameters']['RelativeTolerance'] =    1e-5

    # OpenMP evaluator.
    options['Model']['Evaluators']['Device_0']['Library'] = 'OpenMP'
    options['Model']['Evaluators']['Device_0']['Name']    = 'openmp'
    options['Model']['Evaluators']['Device_0']['Groups']  = ['DiurnalKinetics_C1', 'DiurnalKinetics_C2']
    options['Model']['Evaluators']['Device_0']['Parameters']['numThreads'] = 0

    # ILU options for Ncpu = 1: k = 1, rho = 1.0, alpha = 1e-5, w = 0.0
    options['LinearSolver']['Preconditioner']['Parameters']['fact: level-of-fill']      =    1
    options['LinearSolver']['Preconditioner']['Parameters']['fact: relax value']        =  0.0
    options['LinearSolver']['Preconditioner']['Parameters']['fact: absolute threshold'] = 1e-5
    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, ['C1(0,0)'])
    
if __name__ == "__main__":
    if len(sys.argv) == 1:
        Nx = 80
        Ny = 80
    elif len(sys.argv) == 3:
        Nx = int(sys.argv[1])
        Ny = int(sys.argv[2])
    else:
        print('Usage: python tutorial_opencs_ode_3_groups.py Nx Ny')
        sys.exit()
        
    inputFilesDirectory = 'tutorial_opencs_ode_3_groups'
    run(Nx = Nx, Ny = Ny, inputFilesDirectory = inputFilesDirectory)