#!/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, csUserDefinedKernelFunction 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 dk_kernel_functions = """ #ifndef CS_DIURNAL_KINETICS_UDF_H #define CS_DIURNAL_KINETICS_UDF_H #define __Ny__ %u // Kv is dependent on y-coordinate. static inline real_t Kv(uint32_t y) { const real_t y_0 = %.19f; const real_t y_1 = %.19f; const real_t Kv0 = %.19f; const real_t dy = (y_1 - y_0) / (__Ny__-1); return Kv0 * exp( (y_0 + y * dy) / 5.0 ); } // Kv function used when evaluating derivatives. // It returns the same value but kernels generated in C need the return value to be adouble type. // Therefore, to make the kernel operable in both C and C++ return _constant_(Kv) which returns adouble object. static inline adouble Kv_deriv(uint32_t y) { real_t kv = Kv(y); return _constant_(kv); } #endif """ class DiurnalKineticsKernels_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_BCs = csGroup_t ("BoundaryConditions", 1) # the group for boundary conditions (group 1) self.kernel_C1 = csKernel_t("DiurnalKinetics_C1", 2) # the kernel for C1-component (group 2) self.kernel_C2 = csKernel_t("DiurnalKinetics_C2", 3) # the kernel for C2-component (group 3) self.userDefinedSourceCode = dk_kernel_functions % (self.Ny, self.y0, self.y1, Kv0) 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 csNumber_t(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) apis = csKernel_t.GetSupportedKernelAPIs() for api in apis: udkd_c1 = self.kernel_C1.GetUserDefinedBuildData(api) udkd_c2 = self.kernel_C2.GetUserDefinedBuildData(api) udkd_c1.userDefinedSourceCode = self.userDefinedSourceCode udkd_c2.userDefinedSourceCode = self.userDefinedSourceCode equations = [] kernels = [self.kernel_C1, self.kernel_C2] # Component 2 (C1): for x in range(Nx): for y in range(Ny): if x == 0 or x == Nx-1 or y == 0 or y == Ny-1: equation = csEquation_t(self.group_BCs) 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) else: # GRID_POINT_2D macro takes the csGridPoint_2D_t object from the kernel, set by a call to SetGridPoint(x, y). # It contains two members uint32_t dim_1 and dim_2 representing x, y indexes in two dimensions. fnKv = csUserDefinedKernelFunction(Kv(y), "Kv(GRID_POINT_2D.dim_2)", "Kv_deriv(GRID_POINT_2D.dim_2)") equation_C1 = csEquation_t(self.kernel_C1) equation_C1[ C1(x,y) ] = V * dC1_dx(x,y) + \ Kh * d2C1_dx2(x,y) + \ fnKv * (0.2 * dC1_dy(x,y) + d2C1_dy2(x,y)) + \ R1(C1(x,y), C2(x,y), time) equation_C1.SetGridPoint(x, y); self.kernel_C1.AddEquation(equation_C1); # Component 2 (C2): for x in range(Nx): for y in range(Ny): if x == 0 or x == Nx-1 or y == 0 or y == Ny-1: equation = csEquation_t(self.group_BCs) 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) else: # GRID_POINT_2D macro takes the csGridPoint_2D_t object from the kernel, set by a call to SetGridPoint(x, y). # It contains two members uint32_t dim_1 and dim_2 representing x, y indexes in two dimensions. fnKv = csUserDefinedKernelFunction(Kv(y), "Kv(GRID_POINT_2D.dim_2)", "Kv_deriv(GRID_POINT_2D.dim_2)") equation_C2 = csEquation_t(self.kernel_C2) equation_C2[ C2(x,y) ] = V * dC2_dx(x,y) + \ Kh * d2C2_dx2(x,y) + \ fnKv * (0.2 * dC2_dy(x,y) + d2C2_dy2(x,y)) + \ R2(C1(x,y), C2(x,y), time) equation_C2.SetGridPoint(x, y); self.kernel_C2.AddEquation(equation_C2); 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 = DiurnalKineticsKernels_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 # (1) Sequential evaluator (evaluates groups one by one). options['Model']['Evaluators']['Device_0']['Library'] = 'Sequential' options['Model']['Evaluators']['Device_0']['Name'] = 'seq' options['Model']['Evaluators']['Device_0']['Groups'] = ['BoundaryConditions'] # (2a) OpenMP evaluator. options['Model']['Evaluators']['Device_1'] = {} options['Model']['Evaluators']['Device_1']['Parameters'] = {} options['Model']['Evaluators']['Device_1']['Library'] = 'Kernels_OpenMP' options['Model']['Evaluators']['Device_1']['Name'] = 'kopenmp' options['Model']['Evaluators']['Device_1']['Groups'] = [] options['Model']['Evaluators']['Device_1']['Kernels'] = ['DiurnalKinetics_C1', 'DiurnalKinetics_C2'] options['Model']['Evaluators']['Device_1']['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_kernels.py Nx Ny') sys.exit() inputFilesDirectory = 'tutorial_opencs_ode_3_kernels' run(Nx = Nx, Ny = Ny, inputFilesDirectory = inputFilesDirectory)