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

"""
***********************************************************************************
tutorial_opencs_dae_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 IDAS idasBruss_kry_bbd_p example.
The PDE system is a two-species time-dependent PDE known as
Brusselator PDE and models a chemically reacting system::

du/dt = eps1(d2u/dx2  + d2u/dy2) + u^2 v - (B+1)u + A
dv/dt = eps2(d2v/dx2  + d2v/dy2) - u^2 v + Bu

Boundary conditions: Homogenous Neumann.
Initial Conditions::

u(x,y,t0) = u0(x,y) =  1  - 0.5*cos(pi*y/L)
v(x,y,t0) = v0(x,y) = 3.5 - 2.5*cos(pi*x/L)

The PDEs are discretized by central differencing on a uniform (Nx, Ny) grid.
The model is described in:

- R. Serban and A. C. Hindmarsh. CVODES, the sensitivity-enabled ODE solver in SUNDIALS.
In Proceedings of the 5th International Conference on Multibody Systems,
Nonlinear Dynamics and Control, Long Beach, CA, 2005. ASME.
- M. R. Wittman. Testing of PVODE, a Parallel ODE Solver.
Technical Report UCRL-ID-125562, LLNL, August 1996.

The original results are in tutorial_opencs_dae_3.csv file.
"""

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

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)

eq = 0
equations = [None] * self.Nequations
# Component u:
for x in range(Nx):
for y in range(Ny):
if(x == 0):       # Left BC: Neumann BCs
equations[eq] = du_dx(x,y) - u_flux_bc
elif(x == Nx-1):  # Right BC: Neumann BCs
equations[eq] = du_dx(x,y) - u_flux_bc
elif(y == 0):     # Bottom BC: Neumann BCs
equations[eq] = du_dy(x,y) - u_flux_bc
elif(y == Ny-1):  # Top BC: Neumann BCs
equations[eq] = du_dy(x,y) - u_flux_bc
else:
# Interior points
equations[eq] = 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)
eq += 1

# Component v:
for x in range(Nx):
for y in range(Ny):
if(x == 0):       # Left BC: Neumann BCs
equations[eq] = dv_dx(x,y) - v_flux_bc
elif(x == Nx-1):  # Right BC: Neumann BCs
equations[eq] = dv_dx(x,y) - v_flux_bc
elif(y == 0):     # Bottom BC: Neumann BCs
equations[eq] = dv_dy(x,y) - v_flux_bc
elif(y == Ny-1):  # Top BC: Neumann BCs
equations[eq] = dv_dy(x,y) - v_flux_bc
else:
# Interior points
equations[eq] = 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)
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

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 = Brusselator_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.1
options['Solver']['Parameters']['RelativeTolerance'] = 1e-5

# 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. Generate models for parallel simulations on message-passing multi-processors.
# ILU options for Ncpu = 8: k = 1, rho = 1.0, alpha = 1e-1, 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-1
options['LinearSolver']['Preconditioner']['Parameters']['fact: relative threshold'] =  1.0
mb.SimulationOptions = options

# For distributed memory systems a graph partitioner must be specified.
# Available partitioners:
#   1. csGraphPartitioner_Simple (splits the given set of equations into Npe parts with no analysis)
#   2. csGraphPartitioner_2D_Npde (distributes Npde equations on an uniform 2D grid split into Npe quadrants)
#   3. csGraphPartitioner_Metis (partitions the set of equations using the METIS graph partitioner)
#      Two algorithms are avaiable:
#       - PartGraphKway:      'Multilevel k-way partitioning' algorithm
#       - PartGraphRecursive: 'Multilevel recursive bisectioning' algorithm
#   4. User-defined graph partitioners
# The Metis partitioner can additionally balance the specified quantities in all partitions
# using the following balancing constraints:
#  - Ncs:      balance number of compute stack items in equations
#  - Nnz:      balance number of non-zero items in the incidence matrix
#  - Nflops:   balance number of FLOPS required to evaluate equations
#  - Nflops_j: balance number of FLOPS required to evaluate derivatives (Jacobian) matrix
# PartitionSystem arguments:
#  - Npe: unsigned int specifying the number of processing elements (CPUs)
#  - graphPartitioner: csGraphPartitioner_t instance
#  - balancingConstraints: a list of balancing constraints as strings
#  - logPartitionResults: bool (default False); if true a log file is created
#  - unaryOperationsFlops: dictionary [csUnaryFunctions : unsigned int])
#  - binaryOperationsFlops: dictionary [csBinaryFunctions : unsigned int]
Npe = 8
graphPartitioner = createGraphPartitioner_2D_Npde(Nx = model.Nx,
Ny = model.Ny,
Npde = 2,
Npex_Npey_ratio = 2.0)
inputFilesDirectory_mpi = inputFilesDirectory + ('-Npe=%d-2D_Npde' % Npe)
cs_models_mpi = mb.PartitionSystem(Npe, graphPartitioner, logPartitionResults = True)
csModelBuilder_t.ExportModels(cs_models_mpi, inputFilesDirectory_mpi, mb.SimulationOptions)
print("Npe = %d CPUs OpenCS model generated successfully!" % Npe)

# 5. 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.py Nx Ny')
sys.exit()

u_flux_bc = 0.0
v_flux_bc = 0.0
inputFilesDirectory = 'tutorial_opencs_dae_3'
run(Nx = Nx,
Ny = Ny,
u_flux_bc = u_flux_bc,
v_flux_bc = v_flux_bc,
inputFilesDirectory = inputFilesDirectory)
```