#!/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 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 class DiurnalKineticsGroups_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_C1 = csGroup_t("DiurnalKinetics_C1", 1) # the group for C1-component (group 2) self.group_C2 = csGroup_t("DiurnalKinetics_C2", 2) # the group for C2-component (group 3) 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 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) equations = [] kernels = [] # Component 2 (C1): for x in range(Nx): for y in range(Ny): equation = csEquation_t(self.group_C1) 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) # Component 2 (C2): for x in range(Nx): for y in range(Ny): equation = csEquation_t(self.group_C2) 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) 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 = DiurnalKineticsGroups_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 # OpenMP evaluator. options['Model']['Evaluators']['Device_0']['Library'] = 'OpenMP' options['Model']['Evaluators']['Device_0']['Name'] = 'openmp' options['Model']['Evaluators']['Device_0']['Groups'] = ['DiurnalKinetics_C1', 'DiurnalKinetics_C2'] options['Model']['Evaluators']['Device_0']['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_groups.py Nx Ny') sys.exit() inputFilesDirectory = 'tutorial_opencs_ode_3_groups' run(Nx = Nx, Ny = Ny, inputFilesDirectory = inputFilesDirectory)