#!/usr/bin/env python # -*- coding: utf-8 -*- """ *********************************************************************************** tutorial_cv_11.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__ = """ Code verification using the Method of Exact Solutions. The problem is identical to the problem in the tutorial_cv_6. The only difference is that the flow is reversed to test the high resolution scheme for the reversed flow mode. Numerical vs. exact solution plots (Nx = [20, 40, 80]): .. image:: _static/tutorial_cv_11-results.png :width: 800px """ import sys, math, numpy from time import localtime, strftime import matplotlib.pyplot as plt from daetools.pyDAE import * # Standard variable types are defined in variable_types.py from pyUnits import m, g, kg, s, K, mol, kmol, J, um c_t = daeVariableType("c_t", dimless, -1.0e+20, 1.0e+20, 0.0, 1e-07) # Reversed flow: u < 0 u = -1.0 D = 0.002 L = 1.0 dt = 0.3 pi = numpy.pi class modTutorial(daeModel): def __init__(self, Name, Parent = None, Description = ""): daeModel.__init__(self, Name, Parent, Description) self.x = daeDomain("x", self, m, "") self.c = daeVariable("c", c_t, self, "c using high resolution upwind scheme", [self.x]) def DeclareEquations(self): daeModel.DeclareEquations(self) xp = self.x.Points Nx = self.x.NumberOfPoints t = Time() hr = daeHRUpwindSchemeEquation(self.c, self.x, daeHRUpwindSchemeEquation.Phi_Koren, 1e-10, True) c = lambda i: self.c(i) # Convection-diffusion equation for i in range(0, Nx-1): eq = self.CreateEquation("c(%d)" % i, "") eq.Residual = hr.dc_dt(i) + u * hr.dc_dx(i) - D * hr.d2c_dx2(i) eq.CheckUnitsConsistency = False # BCs eq = self.CreateEquation("c(%d)" % 0, "") eq.Residual = c(Nx-1) - 0.0 class simTutorial(daeSimulation): def __init__(self, Nx): daeSimulation.__init__(self) self.m = modTutorial("tutorial_cv_11(%d)" % Nx) self.m.Description = __doc__ self.Nx = Nx def SetUpParametersAndDomains(self): self.m.x.CreateStructuredGrid(self.Nx, 0.0, L) def SetUpVariables(self): Nx = self.m.x.NumberOfPoints xp = self.m.x.Points for i in range(1, Nx): self.m.c.SetInitialCondition(i, numpy.sin(pi*xp[i])) # Setup everything manually and run in a console def simulate(Nx): # Create Log, Solver, DataReporter and Simulation object log = daePythonStdOutLog() daesolver = daeIDAS() datareporter = daeDelegateDataReporter() simulation = simTutorial(Nx) # Do no print progress log.PrintProgress = False daesolver.RelativeTolerance = 1e-7 # Enable reporting of all variables simulation.m.SetReportingOn(True) # Set the time horizon and the reporting interval simulation.ReportingInterval = 0.02 simulation.TimeHorizon = dt # Connect data reporter simName = simulation.m.Name + strftime(" [%d.%m.%Y %H:%M:%S]", localtime()) # 1. TCP/IP tcpipDataReporter = daeTCPIPDataReporter() datareporter.AddDataReporter(tcpipDataReporter) if not tcpipDataReporter.Connect("", simName): sys.exit() # 2. Data dr = daeNoOpDataReporter() datareporter.AddDataReporter(dr) # Initialize the simulation simulation.Initialize(daesolver, datareporter, log) # Save the model report and the runtime model report #simulation.m.SaveModelReport(simulation.m.Name + ".xml") #simulation.m.SaveRuntimeModelReport(simulation.m.Name + "-rt.xml") # Solve at time=0 (initialization) simulation.SolveInitial() # Run simulation.Run() simulation.Finalize() ########################################### # Data # ########################################### results = dr.Process.dictVariables cvar = results[simulation.m.Name + '.c'] c = cvar.Values[-1, :] # 2D array [t,x] return simulation.m.x.Points, c def run(**kwargs): Nxs = numpy.array([20, 40, 80]) n = len(Nxs) hs = L / Nxs c = [] # Run simulations for i,Nx in enumerate(Nxs): nx, c_ = simulate(int(Nx)) # Exact solution: cexact_ = [] for xk in nx: if xk <= L - numpy.abs(u)*dt: ce = numpy.exp(-D*dt) * numpy.sin(pi*(xk-u*dt)) else: ce = 0 cexact_.append(ce) c.append((nx, c_, cexact_)) fontsize = 14 fontsize_legend = 11 fig = plt.figure(figsize=(12,4), facecolor='white') fig.canvas.set_window_title('Tutorial cv_11') ax = plt.subplot(131) plt.figure(1, facecolor='white') plt.plot(c[0][0], c[0][1], 'ro', linewidth=1.0, label='c (Nx=%d)' % Nxs[0]) plt.plot(c[0][0], c[0][2], 'b-', linewidth=1.0, label='c_exact (Nx=%d)' % Nxs[0]) plt.xlabel('x', fontsize=fontsize) plt.ylabel('c', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) #plt.xlim((0.0, 1.0)) #plt.ylim((0.0, 1.0)) ax = plt.subplot(132) plt.figure(1, facecolor='white') plt.plot(c[1][0], c[1][1], 'ro', linewidth=1.0, label='c (Nx=%d)' % Nxs[1]) plt.plot(c[1][0], c[1][2], 'b-', linewidth=1.0, label='c_exact (Nx=%d)' % Nxs[1]) plt.xlabel('x', fontsize=fontsize) plt.ylabel('c', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) #plt.xlim((0.0, 1.0)) #plt.ylim((0.0, 1.0)) ax = plt.subplot(133) plt.figure(1, facecolor='white') plt.plot(c[2][0], c[2][1], 'ro', linewidth=1.0, label='c (Nx=%d)' % Nxs[2]) plt.plot(c[2][0], c[2][2], 'b-', linewidth=1.0, label='c_exact (Nx=%d)' % Nxs[2]) plt.xlabel('x', fontsize=fontsize) plt.ylabel('c', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) #plt.xlim((0.0, 1.0)) #plt.ylim((0.0, 1.0)) plt.tight_layout() plt.show() if __name__ == "__main__": run()