#!/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.5 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. options['Model']['Evaluators']['Device_0']['Library'] = 'Sequential' options['Model']['Evaluators']['Device_0']['Name'] = 'seq' options['Model']['Evaluators']['Device_0']['Groups'] = ['DefaultGroup'] options['Model']['Evaluators']['Device_0']['Parameters'] = {} # 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)