#!/usr/bin/env python # -*- coding: utf-8 -*- """ *********************************************************************************** tutorial9.py DAE Tools: pyDAE 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__ = """ This tutorial introduces the following concepts: - Third party **direct** linear equations solvers Currently there are the following linear equations solvers available: - SuperLU: sequential sparse direct solver defined in pySuperLU module (BSD licence) - SuperLU_MT: multi-threaded sparse direct solver defined in pySuperLU_MT module (BSD licence) - Trilinos Amesos: sequential sparse direct solver defined in pyTrilinos module (GNU Lesser GPL) - IntelPardiso: multi-threaded sparse direct solver defined in pyIntelPardiso module (proprietary) - Pardiso: multi-threaded sparse direct solver defined in pyPardiso module (proprietary) In this example we use the same conduction problem as in the tutorial 1. The temperature plot (at t=100s, x=0.5, y=*): .. image:: _static/tutorial9-results.png :width: 500px """ import sys from time import localtime, strftime from daetools.pyDAE import * # First import desired solver's module: from daetools.solvers.trilinos import pyTrilinos #from daetools.solvers.superlu import pySuperLU #from daetools.solvers.superlu_mt import pySuperLU_MT #from daetools.solvers.intel_pardiso import pyIntelPardiso #from daetools.solvers.pardiso import pyPardiso # Standard variable types are defined in variable_types.py from pyUnits import m, kg, s, K, J, W class modTutorial(daeModel): def __init__(self, Name, Parent = None, Description = ""): daeModel.__init__(self, Name, Parent, Description) self.x = daeDomain("x", self, m, "X axis domain") self.y = daeDomain("y", self, m, "Y axis domain") self.Qb = daeParameter("Q_b", W/(m**2), self, "Heat flux at the bottom edge of the plate") self.Tt = daeParameter("T_t", K, self, "Temperature at the top edge of the plate") self.rho = daeParameter("ρ", kg/(m**3), self, "Density of the plate") self.cp = daeParameter("c_p", J/(kg*K), self, "Specific heat capacity of the plate") self.k = daeParameter("λ_p", W/(m*K), self, "Thermal conductivity of the plate") self.T = daeVariable("T", temperature_t, self, "Temperature of the plate") self.T.DistributeOnDomain(self.x) self.T.DistributeOnDomain(self.y) def DeclareEquations(self): daeModel.DeclareEquations(self) eq = self.CreateEquation("HeatBalance", "Heat balance equation valid on the open x and y domains") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eOpenOpen) eq.Residual = self.rho() * self.cp() * dt(self.T(x,y)) - self.k() * \ (d2(self.T(x,y), self.x) + d2(self.T(x,y), self.y)) eq = self.CreateEquation("BC_bottom", "Neumann boundary conditions at the bottom edge (constant flux)") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eLowerBound) eq.Residual = - self.k() * d(self.T(x,y), self.y, eCFDM) - self.Qb() eq = self.CreateEquation("BC_top", "Dirichlet boundary conditions at the top edge (constant temperature)") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eUpperBound) eq.Residual = self.T(x,y) - self.Tt() eq = self.CreateEquation("BC_left", "Neumann boundary conditions at the left edge (insulated)") x = eq.DistributeOnDomain(self.x, eLowerBound) y = eq.DistributeOnDomain(self.y, eClosedClosed) eq.Residual = d(self.T(x,y), self.x, eCFDM) eq = self.CreateEquation("BC_right", "Neumann boundary conditions at the right edge (insulated)") x = eq.DistributeOnDomain(self.x, eUpperBound) y = eq.DistributeOnDomain(self.y, eClosedClosed) eq.Residual = d(self.T(x,y), self.x, eCFDM) class simTutorial(daeSimulation): def __init__(self): daeSimulation.__init__(self) self.m = modTutorial("tutorial9") self.m.Description = __doc__ def SetUpParametersAndDomains(self): self.m.x.CreateStructuredGrid(20, 0, 0.1) self.m.y.CreateStructuredGrid(20, 0, 0.1) self.m.k.SetValue(401 * W/(m*K)) self.m.cp.SetValue(385 * J/(kg*K)) self.m.rho.SetValue(8960 * kg/(m**3)) self.m.Qb.SetValue(1e6 * W/(m**2)) self.m.Tt.SetValue(300 * K) def SetUpVariables(self): for x in range(1, self.m.x.NumberOfPoints - 1): for y in range(1, self.m.y.NumberOfPoints - 1): self.m.T.SetInitialCondition(x, y, 300 * K) def run(**kwargs): simulation = simTutorial() # The default linear solver is Sundials dense sequential solver (LU decomposition). # The following 3rd party direct linear solvers are supported: # - Pardiso (multi-threaded - OMP) # - IntelPardiso (multi-threaded - OMP) # - SuperLU (sequential) # - SuperLU_MT (multi-threaded - pthreads, OMP) # - Trilinos Amesos (sequential): Klu, SuperLU, Lapack, Umfpack # If you are using Pardiso or IntelPardiso solvers you have to export their bin directories, # using, for instance, LD_LIBRARY_PATH shell variable (for more details see their documentation). # If you are using OpenMP capable solvers you should set the number of threads # (typically to the number of cores), for instance: # export OMP_NUM_THREADS=4 # or if using IntelPardiso solver: # export MKL_NUM_THREADS=24 # You can place the above command at the end of $HOME/.bashrc (or type it in shell, before simulation). # Import desired solver module (uncomment it from below) and set it using SetLASolver function: #print("Supported Trilinos 3rd party LA solvers:", str(pyTrilinos.daeTrilinosSupportedSolvers())) lasolver = pyTrilinos.daeCreateTrilinosSolver("Amesos_Klu", "") #lasolver = pyTrilinos.daeCreateTrilinosSolver("Amesos_Lapack", "") #lasolver = pyTrilinos.daeCreateTrilinosSolver("Amesos_Umfpack", "") #lasolver = pyIntelPardiso.daeCreateIntelPardisoSolver() #lasolver = pyPardiso.daeCreatePardisoSolver() """ # Get Pardiso/IntelPardiso parameters (iparm[64] list of integers) iparm = lasolver.get_iparm() iparm_def = list(iparm) # Save it for comparison # Change some options # Achtung, Achtung!! # The meaning of items in iparm[64] is NOT identical in Pardiso and IntelPardiso solvers!! iparm[ 7] = 2 # Max. number of iterative refinement steps (common for both Pardiso and IntelPardiso) #iparm[27] = 1 # in Pardiso it means: use METIS parallel reordering # in IntelPardiso it means: use single precision (do not change it!) # Set them back lasolver.set_iparm(iparm) iparm = lasolver.get_iparm() # Print the side by side comparison print('iparm default new') for i in range(64): print 'iparm[%2d] %7d %3d' % (i, iparm_def[i], iparm[i]) """ return daeActivity.simulate(simulation, reportingInterval = 10, timeHorizon = 1000, lasolver = lasolver, **kwargs) if __name__ == "__main__": guiRun = False if (len(sys.argv) > 1 and sys.argv[1] == 'console') else True run(guiRun = guiRun)