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

"""
***********************************************************************************
tutorial_dealii_8.py
DAE Tools: pyDAE 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__ = """
In this example a small parallel-plate reactor with an active surface is modelled.

This problem and its solution in `COMSOL Multiphysics <https://www.comsol.com>`_ software
is described in the Application Gallery: `Transport and Adsorption (id=5)

The transport in the bulk of the reactor is described by a convection-diffusion equation:

.. code-block:: none

dc/dt - D*nabla^2(c) + div(uc) = 0 in Omega

The material balance for the surface, including surface diffusion and the reaction rate is:

.. code-block:: none

dc_s/dt - Ds*nabla^2(c_s) = -k_ads * c * (Gamma_s - c_s) + k_des * c_s in Omega_s

For the bulk, the boundary condition at the active surface couples the rate of the reaction
at the surface with the flux of the reacting species and the concentration of the adsorbed
species and bulk species:

.. code-block:: none

n⋅(-D*nabla(c) + c*u) = -k_ads*c*(Gamma_s - c_s) + k_des*c_s

The boundary conditions for the surface species are insulating conditions.

The problem is modelled using two coupled Finite Element systems:
2D for bulk flow and 1D for the active surface.
The linear interpolation is used to determine the bulk flow and active surface concentrations
at the interface.

The mesh is rectangular with the refined elements close to the left/right ends:

.. image:: _static/parallel_plate_reactor.png
:height: 400 px

The cs plot at t = 2s:

.. image:: _static/tutorial_dealii_8-results.png
:height: 400 px

The c plot at t = 2s:

.. image:: _static/tutorial_dealii_8-results2.png
:height: 400 px

Animation:

.. image:: _static/tutorial_dealii_8-animation.gif
:height: 400 px
"""

import os, sys, numpy, scipy.interpolate, operator, tempfile
from time import localtime, strftime
from daetools.pyDAE import *
from daetools.solvers.deal_II import *
from daetools.solvers.superlu import pySuperLU

# Standard variable types are defined in variable_types.py
from pyUnits import m, kg, s, K, Pa, mol, J, W

# Inputs:
c0      = 1000  # [mol/m^3]     Initial concentration
k_ads   = 1e-6  # [m^3/(mol*s)] Forward rate constant
k_des   = 1e-9  # [1/s]	        Backward rate constant
Gamma_s = 1000  # [mol/m^2]     Active site concentration
Ds      = 1e-11 # [m^2/s]       Surface diffusivity
D       = 1e-9  # [m^2/s]       Gas diffusivity
v_max   = 1e-3  # [m/s]         Maximum velocity
delta   = 1e-4  # [m]           Channel width

# Boundary IDs (reactor)
left_edge      = 0
top_edge       = 1
right_edge     = 2 # excluding the active surface
bottom_edge    = 3
active_surface = 5

# Boundary IDs (active surface)
bottom_point = 0
top_point    = 1

dx = 1e-10
dy = 1e-10
def points_equal(x1, y1, x2, y2):
return (x1 < x2+dx and x1 > x2-dx and y1 < y2+dy and y1 > y2-dy)

def coord_equal(x1, x2):
#print(x1, x2, (x1 < x2+dx and x1 > x2-dx))
return (x1 < x2+dx and x1 > x2-dx)

class VelocityFunction_2D(Function_2D):
def __init__(self, n_components = 1):
Function_2D.__init__(self, n_components)
self.u = Tensor_1_2D()

def gradient(self, point, component = 0):
self.u[0] = 0.0
self.u[1] = v_max * (1 - ((point.x-0.5*delta)/(0.5*delta))**2)
return self.u

return [self.value(point, c) for c in range(self.n_components)]

class c_interpolation(object):
"""
Returns linear interpolation of the concentration based on y coordinates.
"""
def __init__(self, c_indexes):
c_i = sorted(c_indexes.items(), key = operator.itemgetter(1))
n = len(c_i)

self.ci = numpy.empty((n), dtype=numpy.int_)    # array of indexes
self.y  = numpy.empty((n), dtype=numpy.float_)  # array of y coordinates
self.c  = numpy.empty((n), dtype=object)        # array of adouble objects

for i, (ci, (y, c)) in enumerate(c_i):
self.ci[i] = ci
self.y[i]  = y
self.c[i]  = c

def get_c(self, y):
if y < self.y[0] or y > self.y[-1]:
print('y = %f out of bounds' % y)
elif y == self.y[0]:
c = self.c[0]
elif y == self.y[-1]:
c = self.c[-1]
else:
i0 = numpy.argwhere(self.y <= y)
if len(i0):
index = i0[-1][0]
y0 = self.y[index]
c0 = self.c[index]
else:

i1 = numpy.argwhere(self.y >= y)
if len(i1):
index = i1[0][0]
y1 = self.y[index]
c1 = self.c[index]
else:

c = (c0*(y1-y) + c1*(y-y0)) / (y1-y0)

return c

def __init__(self, c_indexes, n_components = 1):
c_interpolation.__init__(self, c_indexes)

def value(self, point, component = 0):
return self.get_c(point.y)

def vector_value(self, point):
return [self.value(point, c) for c in range(self.n_components)]

def __init__(self, c_indexes, n_components = 1):
c_interpolation.__init__(self, c_indexes)

def value(self, point, component = 0):
return self.get_c(point.x)

def vector_value(self, point):
return [self.value(point, c) for c in range(self.n_components)]

class modTutorial(daeModel):
def __init__(self, Name, Parent = None, Description = ""):
daeModel.__init__(self, Name, Parent, Description)

self.FE_c_init()
self.FE_cs_init()

def FE_c_init(self):
# FE system/model for the bulk gas concentration.
dofs = [dealiiFiniteElementDOF_2D(name='c',
description='Bulk gas concentration',
fe = FE_Q_2D(1),
multiplicity=1)]

meshes_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'meshes')
mesh_file  = os.path.join(meshes_dir, 'parallel_plate_reactor.msh')

self.fe_system_c = dealiiFiniteElementSystem_2D(meshFilename    = mesh_file,     # path to mesh
dofs            = dofs)          # degrees of freedom

self.fe_model_c = daeFiniteElementModel('BulkFlow', self, 'Transient convection-diffusion-reaction equation', self.fe_system_c)

self.c = self.fe_model_c.Variables[0]

def FE_cs_init(self):
# FE system/model for the active surface.
dofs = [dealiiFiniteElementDOF_1D(name         = 'cs',
description  = 'Surface concentration',
fe           = FE_Q_1D(1),
multiplicity = 1)]

meshes_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'meshes')
mesh_file  = os.path.join(meshes_dir, 'active_surface.msh')

self.fe_system_cs = dealiiFiniteElementSystem_1D(meshFilename    = mesh_file,     # path to mesh
dofs            = dofs)          # degrees of freedom

self.fe_model_cs = daeFiniteElementModel('ActiveSurface', self, 'Transient diffusion-reaction equation', self.fe_system_cs)

self.cs = self.fe_model_cs.Variables[0]

def DeclareEquations(self):
daeModel.DeclareEquations(self)

# Populate dictionaries:
# c_indexes = {overallIndex : (y-coord,  c(overallIndex))}
# c_indexes = {overallIndex : (y-coord, cs(overallIndex))}
# They will be used to get the concentration c/cs by linear interpolation.
self.c_indexes = {}
self.cs_indexes = {}

c_sp  = list(self.fe_system_c.GetDOFSupportPoints())
cs_sp = list(self.fe_system_cs.GetDOFSupportPoints())
for ci, pcs in enumerate(cs_sp):
self.cs_indexes[ci] = (pcs.x, self.cs(ci))
for ci, pc in enumerate(c_sp):
if pc.x >= delta - dx and pc.y >= 0.0 and pc.y <= delta + dy:
self.c_indexes[ci] = (pc.y, self.c(ci))

# Create weak formulations
self.WeakForm_c()
self.WeakForm_cs()

def WeakForm_c(self):
# Create some auxiliary objects for readability
phi_i  =  phi_2D('c', fe_i, fe_q)
phi_j  =  phi_2D('c', fe_j, fe_q)
dphi_i = dphi_2D('c', fe_i, fe_q)
dphi_j = dphi_2D('c', fe_j, fe_q)
normal = normal_2D(fe_q)
xyz    = xyz_2D(fe_q)
JxW    = JxW_2D(fe_q)
c_dof  = dof_approximation_2D('c', fe_q)

dirichletBC = {}

# Function<dim> wrapper

self.cs_fun = cs_interpolation_2D(self.cs_indexes)

# FE weak form terms
c_accumulation = (phi_i * phi_j) * JxW
c_diffusion    = (dphi_i * dphi_j) * D * JxW
c_convection   = phi_i * (u_grad * dphi_j) * JxW
c_source       = 0.0 * JxW

c_faceFluxes = {}
c_flux = -k_ads * c_dof * (Gamma_s - cs_fun) + k_des * cs_fun
c_faceFluxes[active_surface] = phi_i * c_flux * JxW
#c_faceFluxes[top_edge] = phi_i * (normal * u_grad) * c_dof * JxW

weakForm = dealiiFiniteElementWeakForm_2D(Aij = c_diffusion + c_convection,
Mij = c_accumulation,
Fi  = c_source,
boundaryFaceFi = c_faceFluxes,
functionsDirichletBC = dirichletBC)

self.fe_system_c.WeakForm = weakForm

def WeakForm_cs(self):
# Create some auxiliary objects for readability
phi_i  =  phi_1D('cs', fe_i, fe_q)
phi_j  =  phi_1D('cs', fe_j, fe_q)
dphi_i = dphi_1D('cs', fe_i, fe_q)
dphi_j = dphi_1D('cs', fe_j, fe_q)
xyz    = xyz_1D(fe_q)
JxW    = JxW_1D(fe_q)
cs_dof = dof_approximation_1D('cs', fe_q)

self.c_fun = c_interpolation_1D(self.c_indexes)

# FE weak form terms
cs_accumulation = (phi_i * phi_j) * JxW
cs_diffusion    = (dphi_i * dphi_j) * Ds * JxW
cs_flux         = k_ads * c_fun * (Gamma_s - cs_dof) - k_des * cs_dof
cs_source       = phi_i * cs_flux * JxW

weakForm = dealiiFiniteElementWeakForm_1D(Aij = cs_diffusion,
Mij = cs_accumulation,
Fi  = cs_source)

self.fe_system_cs.WeakForm = weakForm

class simTutorial(daeSimulation):
def __init__(self):
daeSimulation.__init__(self)
self.m = modTutorial("tutorial_dealii_8")
self.m.Description = __doc__
self.m.fe_model_c.Description  = __doc__
self.m.fe_model_cs.Description = __doc__

def SetUpParametersAndDomains(self):
pass

def SetUpVariables(self):
setFEInitialConditions(self.m.fe_model_c,  self.m.fe_system_c,  'c',  c0)
setFEInitialConditions(self.m.fe_model_cs, self.m.fe_system_cs, 'cs', 0.0)

def run(**kwargs):
guiRun = kwargs.get('guiRun', False)

simulation = simTutorial()

# Create SuperLU LA solver
lasolver = pySuperLU.daeCreateSuperLUSolver()

# Create and setup two data reporters:
datareporter = daeDelegateDataReporter()
simName = simulation.m.Name + strftime(" [%d.%m.%Y %H:%M:%S]", localtime())
if guiRun:
results_folder_c  = tempfile.mkdtemp(suffix = '-results_c',  prefix = 'tutorial_deal_II_8-')
results_folder_cs = tempfile.mkdtemp(suffix = '-results_cs', prefix = 'tutorial_deal_II_8-')
daeQtMessage("deal.II", "The simulation results will be located in: %s and %s" % (results_folder_c, results_folder_cs))
else:
results_folder_c  = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'tutorial_deal_II_8-results-c')
results_folder_cs = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'tutorial_deal_II_8-results-cs')
print("The simulation results will be located in: %s and %s" % (results_folder_c, results_folder_cs))

# 1. deal.II for c
feDataReporter_c = simulation.m.fe_system_c.CreateDataReporter()
if not feDataReporter_c.Connect(results_folder_c, simName):
sys.exit()

# 2. deal.II for cs
feDataReporter_cs = simulation.m.fe_system_cs.CreateDataReporter()
if not feDataReporter_cs.Connect(results_folder_cs, simName):
sys.exit()

# 3. TCP/IP
tcpipDataReporter = daeTCPIPDataReporter()