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

"""
***********************************************************************************
tutorial_opencs_dae_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__ = """
Reimplementation of DAE Tools tutorial1.py example.
A simple heat conduction problem: conduction through a very thin, rectangular copper plate::

rho * cp * dT(x,y)/dt = k * [d2T(x,y)/dx2 + d2T(x,y)/dy2];  x in (0, Lx), y in (0, Ly)

Two-dimensional Cartesian grid (x,y) of 20 x 20 elements.
The original results are in tutorial_opencs_dae_2.csv file.
"""

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

rho = 8960 # density, kg/m^3
cp  =  385 # specific heat capacity, J/(kg.K)
k   =  401 # thermal conductivity, W/(m.K)
Qb  =  1E5 # flux at the bottom edge, W/m^2
Tt  =  300 # T at the top edge, K

class HeatConduction_2D:
def __init__(self, Nx, Ny):
self.Nx = Nx
self.Ny = Ny

self.Lx = 0.1 # m
self.Ly = 0.1 # m

self.dx = self.Lx / (Nx-1)
self.dy = self.Ly / (Ny-1)

self.Nequations = Nx*Ny

def GetInitialConditions(self):
y0 = [300.0] * self.Nequations
return y0

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

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
T_values = y
T_derivs = dydt
dx = self.dx
dy = self.dy
Nx = self.Nx
Ny = self.Ny

def T(x, y):
index = self.GetIndex(x,y)
return T_values[index]

def dT_dt(x, y):
index = self.GetIndex(x,y)
return T_derivs[index]

# First order partial derivative per x.
def dT_dx(x, y):
if (x == 0): # left
T0 = T(0, y)
T1 = T(1, y)
T2 = T(2, y)
return (-3*T0 + 4*T1 - T2) / (2*dx)
elif (x == Nx-1): # right
Tn  = T(Nx-1,   y)
Tn1 = T(Nx-1-1, y)
Tn2 = T(Nx-1-2, y)
return (3*Tn - 4*Tn1 + Tn2) / (2*dx)
else:
T1 = T(x+1, y)
T2 = T(x-1, y)
return (T1 - T2) / (2*dx)

# First order partial derivative per y.
def dT_dy(x, y):
if (y == 0): # bottom
T0 = T(x, 0)
T1 = T(x, 1)
T2 = T(x, 2)
return (-3*T0 + 4*T1 - T2) / (2*dy)
elif (y == Ny-1): # top
Tn  = T(x, Ny-1  )
Tn1 = T(x, Ny-1-1)
Tn2 = T(x, Ny-1-2)
return (3*Tn - 4*Tn1 + Tn2) / (2*dy)
else:
Ti1 = T(x, y+1)
Ti2 = T(x, y-1)
return (Ti1 - Ti2) / (2*dy)

# Second order partial derivative per x.
def d2T_dx2(x, y):
# This function is typically called only for interior points.
if (x == 0 or x == Nx-1):
raise RuntimeError("d2T_dx2 called for boundary x point")

Ti1 = T(x+1, y)
Ti  = T(x,   y)
Ti2 = T(x-1, y)
return (Ti1 - 2*Ti + Ti2) / (dx*dx)

# Second order partial derivative per y.
def d2T_dy2(x, y):
# This function is typically called only for interior points.
if (y == 0 or y == Ny-1):
raise RuntimeError("d2T_dy2 called for boundary y point")

Ti1 = T(x, y+1)
Ti  = T(x,   y)
Ti2 = T(x, y-1)
return (Ti1 - 2*Ti + Ti2) / (dy*dy)

eq = 0
equations = [None] * self.Nequations
for x in range(Nx):
for y in range(Ny):
if (x == 0):                                # Left BC: zero flux
equations[eq] = dT_dx(x,y)
elif (x == Nx-1):                           # Right BC: zero flux
equations[eq] = dT_dx(x,y)
elif (x > 0 and x < Nx-1 and y == 0):       # Bottom BC: prescribed flux
equations[eq] = -k * dT_dy(x,y) - Qb
elif (x > 0 and x < Nx-1 and y == Ny-1):     # Top BC: prescribed flux
equations[eq] = T(x,y) - Tt
else:                                       # Inner region: diffusion equation
equations[eq] = rho * cp * dT_dt(x,y) - k * (d2T_dx2(x,y) + d2T_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', 20)
Ny = kwargs.get('Ny', 20)

# Instantiate the model being simulated.
model = HeatConduction_2D(Nx, Ny)

# 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-10)

# 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
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, ['T(0,0)'])

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

inputFilesDirectory = 'tutorial_opencs_dae_2'
run(Nx = Nx, Ny = Ny, inputFilesDirectory = inputFilesDirectory)
```