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

"""
***********************************************************************************
tutorial_opencs_ode_2.py
DAE Tools: pyOpenCS module, www.daetools.com
***********************************************************************************
DAE Tools is free software; you can redistribute it and/or modify it under the
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__ = """
The problem is simple advection-diffusion in 2-D::

du/dt = d2u/dx2 + 0.5 du/dx + d2u/dy2

on the rectangle::

0 <= x <= 2
0 <= y <= 1

and simulated for 1 s.
Homogeneous Dirichlet boundary conditions are imposed, with the initial conditions::

u(x,y,t=0) = x(2-x)y(1-y)exp(5xy)

The PDE is discretized on a uniform Nx+2 by Ny+2 grid with central differencing.
The boundary points are eliminated leaving an ODE system of size Nx*Ny.
The original results are in tutorial_opencs_ode_2.csv file.
"""

import os, sys, json, itertools, numpy
from daetools.solvers.opencs import csModelBuilder_t, csNumber_t, csSimulate
from daetools.examples.tutorial_opencs_aux import compareResults

def __init__(self, Nx, Ny, u_bc):
#In the CVode example cvsAdvDiff_bnd.c they only modelled interior points,
#  excluded the boundaries from the ODE system, and assumed homogenous Dirichlet BCs (0.0).
#There, they divided the 2D domain into (Nx+1) by (Ny+1) points and
#  the points at x=0, x=Lx, y=0 and y=Ly are not used in the model.
#Thus, x domain starts at x=1*dx, y domain starts at x=1*dy.
self.Nx   = Nx
self.Ny   = Ny
self.u_bc = u_bc

self.x0 = 0.0
self.x1 = 2.0
self.y0 = 0.0
self.y1 = 1.0
self.dx = (self.x1-self.x0) / (self.Nx+2-1)
self.dy = (self.y1-self.y0) / (self.Ny+2-1)

self.Nequations = self.Nx*self.Ny

def GetInitialConditions(self):
u0 = [0.0] * self.Nequations

x0 = self.x0
x1 = self.x1
y0 = self.y0
y1 = self.y1
dx = self.dx
dy = self.dy
for ix in range(self.Nx):
for iy in range(self.Ny):
index = self.GetIndex(ix,iy)
x = (ix+1) * dx
y = (iy+1) * dy
u0[index] = x*(x1 - x)*y*(y1 - y)*numpy.exp(5*x*y)
return u0

def GetVariableNames(self):
return ['u(%d,%d)'%(x,y) for x,y in itertools.product(range(self.Nx), range(self.Ny))]

def CreateEquations(self, y):
# y is a list of csNumber_t objects representing model variables
u_values = y
dx = self.dx
dy = self.dy
Nx = self.Nx
Ny = self.Ny

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

# First order partial derivative per x.
def du_dx(x, y):
# If called for x == 0 or x == Nx-1 use the boundary value (u_bc = 0.0 in this example).
ui1 = (self.u_bc if x == Nx-1 else u(x+1, y))
ui2 = (self.u_bc if x == 0    else u(x-1, y))
return (ui1 - ui2) / (2*dx)

# First order partial derivative per y (not used in this example).
def du_dy(x, y):
# If called for y == 0 or y == Ny-1 use the boundary value (u_bc = 0.0 in this example).
ui1 = (self.u_bc if y == Ny-1 else u(x, y+1))
ui2 = (self.u_bc if y == 0    else u(x, y-1))
return (ui1 - ui2) / (2*dy)

# Second order partial derivative per x.
def d2u_dx2(x, y):
# If called for x == 0 or x == Nx-1 use the boundary value (u_bc = 0.0 in this example).
ui1 = (self.u_bc if x == Nx-1 else u(x+1, y))
ui  =                              u(x,   y)
ui2 = (self.u_bc if x == 0    else u(x-1, y))
return (ui1 - 2*ui + ui2) / (dx*dx)

# Second order partial derivative per y.
def d2u_dy2(x, y):
# If called for y == 0 or y == Ny-1 use the boundary value (u_bc = 0.0 in this example).
ui1 = (self.u_bc if y == Ny-1 else u(x, y+1))
ui  =                              u(x,   y)
ui2 = (self.u_bc if y == 0    else u(x, y-1))
return (ui1 - 2*ui + ui2) / (dy*dy)

eq = 0
equations = [None] * self.Nequations
for x in range(Nx):
for y in range(Ny):
equations[eq] = d2u_dx2(x,y) + 0.5 * du_dx(x,y) + d2u_dy2(x,y)
eq += 1

return equations

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',   10)
Ny   = kwargs.get('Ny',   5)
u_bc = kwargs.get('u_bc', 0.0)

# Instantiate the model being simulated.

# 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-6, defaultVariableName = 'u')

# 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)
# 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']                 =  1.0
options['Simulation']['ReportingInterval']           =  0.1
options['Solver']['Parameters']['RelativeTolerance'] = 1e-5
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("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(9,4)'])

if __name__ == "__main__":
if len(sys.argv) == 1:
Nx = 10
Ny = 5
elif len(sys.argv) == 3:
Nx = int(sys.argv[1])
Ny = int(sys.argv[2])
else:
print('Usage: python tutorial_opencs_ode_2.py Nx Ny')
sys.exit()

u_bc = 0.0
inputFilesDirectory = 'tutorial_opencs_ode_2'
run(Nx = Nx, Ny = Ny, u_bc = u_bc, inputFilesDirectory = inputFilesDirectory)
```