#!/usr/bin/env python # -*- coding: utf-8 -*- """ *********************************************************************************** dae_example_4_cv_plots.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 Manufactured Solutions. The problem in this tutorial is the transient convection-diffusion equation distributed on a rectangular 2D domain with u and v components of velocity (Burgers equations). The problem is identical to DAE Tools Code Verification tutorial_cv_4.py example (also included here for comparison purposes). Testing procedure: 1. Run dae_example_4_cv It will produce dae_example_4_cv-sequential directory with 10x8, 20x16, 40x32, 80x64 subdirectories with results for the corresponding meshes. 2. Run python dae_example_4_cv_plots.py opencs|comparison. (a) Using the opencs option, the script will load results from dae_example_4_cv-sequential and plot the normalised global errors and the order of accuracy plots for OpenCS results. (b) Using comparison option, the script will first run daetools simulations for the 10x8, 20x16, 40x32 and 80x64 meshes, load OpenCS data from the dae_example_4_cv-sequential directory, print the comparison and plot the DAE Tools normalised global errors and the order of accuracy plots. """ import sys, numpy from time import localtime, strftime import matplotlib.pyplot as plt from daetools.pyDAE import * from daetools.solvers.superlu import pySuperLU import pandas no_t = daeVariableType("no_t", dimless, -1.0e+20, 1.0e+20, 0.0, 1e-6) 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.u0 = 1.0 self.v0 = 1.0 self.w0 = 0.1 self.ni = 0.7 self.eps = 0.001 self.u = daeVariable("u", no_t, self, "", [self.x, self.y]) self.v = daeVariable("v", no_t, self, "", [self.x, self.y]) self.um = daeVariable("um", no_t, self, "", [self.x, self.y]) self.vm = daeVariable("vm", no_t, self, "", [self.x, self.y]) def DeclareEquations(self): daeModel.DeclareEquations(self) # Create some auxiliary functions to make equations more readable # DAE Tools can do symbolic differentiation whereas OpenCS cannot. u0 = self.u0 v0 = self.v0 ni = self.ni w0 = self.w0 eps = self.eps t = Time() u = lambda x,y: self.u(x,y) v = lambda x,y: self.v(x,y) du_dt = lambda x,y: dt(u(x,y)) duu_dx = lambda x,y: d(u(x,y)*u(x,y), self.x) duv_dy = lambda x,y: d(u(x,y)*v(x,y), self.y) d2u_dx2 = lambda x,y: d2(u(x,y), self.x) d2u_dy2 = lambda x,y: d2(u(x,y), self.y) dv_dt = lambda x,y: dt(v(x,y)) dvu_dx = lambda x,y: d(v(x,y)*u(x,y), self.x) dvv_dy = lambda x,y: d(v(x,y)*v(x,y), self.y) d2v_dx2 = lambda x,y: d2(v(x,y), self.x) d2v_dy2 = lambda x,y: d2(v(x,y), self.y) um = lambda x,y: u0 * (numpy.sin(x()**2 + y()**2 + w0*t) + eps) dum_dt = lambda x,y: dt(um(x,y)) dumum_dx = lambda x,y: d(um(x,y)*um(x,y), self.x) dumvm_dy = lambda x,y: d(um(x,y)*vm(x,y), self.y) d2um_dx2 = lambda x,y: d2(um(x,y), self.x) d2um_dy2 = lambda x,y: d2(um(x,y), self.y) vm = lambda x,y: v0 * (numpy.cos(x()**2 + y()**2 + w0*t) + eps) dvm_dt = lambda x,y: dt(vm(x,y)) dvmum_dx = lambda x,y: d(vm(x,y)*um(x,y), self.x) dvmvm_dy = lambda x,y: d(vm(x,y)*vm(x,y), self.y) d2vm_dx2 = lambda x,y: d2(vm(x,y), self.x) d2vm_dy2 = lambda x,y: d2(vm(x,y), self.y) Su = lambda x,y: dum_dt(x,y) + (dumum_dx(x,y) + dumvm_dy(x,y)) - ni * (d2um_dx2(x,y) + d2um_dy2(x,y)) Sv = lambda x,y: dvm_dt(x,y) + (dvmum_dx(x,y) + dvmvm_dy(x,y)) - ni * (d2vm_dx2(x,y) + d2vm_dy2(x,y)) # Numerical solution eq = self.CreateEquation("u", "u velocity component") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eOpenOpen) eq.Residual = du_dt(x,y) + (duu_dx(x,y) + duv_dy(x,y)) - ni * (d2u_dx2(x,y) + d2u_dy2(x,y)) - Su(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("u(,0)", "BCs u component") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eLowerBound) eq.Residual = u(x,y) - um(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("u(,1)", "BCs u component") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eUpperBound) eq.Residual = u(x,y) - um(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("u(0,)", "BCs u component") x = eq.DistributeOnDomain(self.x, eLowerBound) y = eq.DistributeOnDomain(self.y, eClosedClosed) eq.Residual = u(x,y) - um(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("u(1,)", "BCs u component") x = eq.DistributeOnDomain(self.x, eUpperBound) y = eq.DistributeOnDomain(self.y, eClosedClosed) eq.Residual = u(x,y) - um(x,y) eq.CheckUnitsConsistency = False # v component eq = self.CreateEquation("v", "v velocity component") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eOpenOpen) eq.Residual = dv_dt(x,y) + (dvu_dx(x,y) + dvv_dy(x,y)) - ni * (d2v_dx2(x,y) + d2v_dy2(x,y)) - Sv(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("v(,0)", "BCs v component") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eLowerBound) eq.Residual = v(x,y) - vm(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("v(,1)", "BCs v component") x = eq.DistributeOnDomain(self.x, eOpenOpen) y = eq.DistributeOnDomain(self.y, eUpperBound) eq.Residual = v(x,y) - vm(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("v(0,)", "BCs v component") x = eq.DistributeOnDomain(self.x, eLowerBound) y = eq.DistributeOnDomain(self.y, eClosedClosed) eq.Residual = v(x,y) - vm(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("v(1,)", "BCs v component") x = eq.DistributeOnDomain(self.x, eUpperBound) y = eq.DistributeOnDomain(self.y, eClosedClosed) eq.Residual = v(x,y) - vm(x,y) eq.CheckUnitsConsistency = False # Manufactured solution eq = self.CreateEquation("um", "u component manufactured solution") x = eq.DistributeOnDomain(self.x, eClosedClosed) y = eq.DistributeOnDomain(self.y, eClosedClosed) eq.Residual = self.um(x,y) - um(x,y) eq.CheckUnitsConsistency = False eq = self.CreateEquation("vm", "v component manufactured solution") x = eq.DistributeOnDomain(self.x, eClosedClosed) y = eq.DistributeOnDomain(self.y, eClosedClosed) eq.Residual = self.vm(x,y) - vm(x,y) eq.CheckUnitsConsistency = False class simTutorial(daeSimulation): def __init__(self, Nx, Ny): daeSimulation.__init__(self) self.m = modTutorial("tutorial_cv_4(%dx%d)" % (Nx,Ny)) self.m.Description = __doc__ self.Nx = Nx self.Ny = Ny def SetUpParametersAndDomains(self): self.m.x.CreateStructuredGrid(self.Nx, -0.1, 0.7) self.m.y.CreateStructuredGrid(self.Ny, 0.2, 0.8) def SetUpVariables(self): Nx = self.m.x.NumberOfPoints Ny = self.m.y.NumberOfPoints xp = self.m.x.Points yp = self.m.y.Points u0 = self.m.u0 v0 = self.m.v0 eps = self.m.eps um0 = lambda x,y: u0 * (numpy.sin(x**2 + y**2) + eps) vm0 = lambda x,y: v0 * (numpy.cos(x**2 + y**2) + eps) for x in range(1, Nx-1): for y in range(1, Ny-1): self.m.u.SetInitialCondition(x,y, um0(xp[x], yp[y])) self.m.v.SetInitialCondition(x,y, vm0(xp[x], yp[y])) # Setup everything manually and run in a console def simulate(Nx, Ny, **kwargs): # Create Log, Solver, DataReporter and Simulation object log = daePythonStdOutLog() daesolver = daeIDAS() datareporter = daeDelegateDataReporter() simulation = simTutorial(Nx, Ny) # Do no print progress log.PrintProgress = False lasolver = pySuperLU.daeCreateSuperLUSolver() # Enable reporting of time derivatives for all reported variables simulation.ReportTimeDerivatives = True # 1. TCP/IP #tcpipDataReporter = daeTCPIPDataReporter() #datareporter.AddDataReporter(tcpipDataReporter) #simName = simulation.m.Name + strftime(" [%d.%m.%Y %H:%M:%S]", localtime()) #if not tcpipDataReporter.Connect("", simName): # sys.exit() # 2. Data dr = daeNoOpDataReporter() datareporter.AddDataReporter(dr) daeActivity.simulate(simulation, reportingInterval = 2.0, timeHorizon = 90.0, saveRuntimeModelReport = True, log = log, lasolver = lasolver, datareporter = datareporter, **kwargs) ########################################### # Data # ########################################### results = dr.Process.dictVariables uvar = results[simulation.m.Name + '.u'] vvar = results[simulation.m.Name + '.v'] umvar = results[simulation.m.Name + '.um'] vmvar = results[simulation.m.Name + '.vm'] times = uvar.TimeValues u = uvar.Values[-1, :, :] # 3D array [t,x,y] v = vvar.Values[-1, :, :] # 3D array [t,x,y] um = umvar.Values[-1, :, :] # 3D array [t,x,y] vm = vmvar.Values[-1, :, :] # 3D array [t,x,y] return times,u,v,um,vm, simulation Nxs_ = [10, 20, 40, 80] Nys_ = [ 8, 16, 32, 64] def loadOpenCSData(resultsDirectory = 'dae_example_4_cv-sequential'): """ Get the data and return the normalised global errors. """ Nxs = numpy.array(Nxs_) Nys = numpy.array(Nys_) Nexperiments = len(Nxs) Lx = 0.8 # Length of the x domain Eu = numpy.zeros(Nexperiments) Ev = numpy.zeros(Nexperiments) u = [0]*Nexperiments v = [0]*Nexperiments um = [0]*Nexperiments vm = [0]*Nexperiments for i in range(Nexperiments): Nx = int(Nxs[i]) Ny = int(Nys[i]) csvFilename = os.path.join(resultsDirectory, '%dx%d' % (Nx, Ny), 'results-0.csv') csv_filepath = os.path.join(os.path.dirname(os.path.abspath(__file__)), csvFilename) df = pandas.read_csv(csv_filepath, sep=';', header=2, skiprows=None, quotechar='"', skipinitialspace=True, dtype=float) #df.info() #print(df.columns) Nx = Nx+1 Ny = Ny+1 print('Loading %dx%d data' % (Nx-1, Ny-1)) u_start = 1 # The first value is Time (skip it) u_end = u_start + Nx*Ny v_start = u_end v_end = v_start + Nx*Ny um_start = v_end um_end = um_start + Nx*Ny vm_start = um_end vm_end = vm_start + Nx*Ny #print(u_start,u_end, v_start,v_end, um_start,um_end, vm_start,vm_end) data = df.values[-1] # take the last time point ul = data[u_start:u_end] vl = data[v_start:v_end] uml = data[um_start:um_end] vml = data[vm_start:vm_end] #if i == 0: # print('%dx%d:' % (Nx, Ny)) # for val in zip(ul - uml, vl - vml): # print(val) u[i] = numpy.array(ul).reshape(Nx,Ny) v[i] = numpy.array(vl).reshape(Nx,Ny) um[i] = numpy.array(uml).reshape(Nx,Ny) vm[i] = numpy.array(vml).reshape(Nx,Ny) Eu[i] = numpy.sqrt((1.0/(Nx*Ny)) * numpy.sum((u[i]-um[i])**2)) Ev[i] = numpy.sqrt((1.0/(Nx*Ny)) * numpy.sum((v[i]-vm[i])**2)) return u, v, um, vm, Eu, Ev def compareResults(u, v, um, vm, u_dt, v_dt, um_dt, vm_dt): Nxs = numpy.array(Nxs_) Nys = numpy.array(Nys_) Nexperiments = len(Nxs) Eu = numpy.zeros(Nexperiments) Ev = numpy.zeros(Nexperiments) Eum = numpy.zeros(Nexperiments) Evm = numpy.zeros(Nexperiments) print('Comparison to daetools data') for i in range(Nexperiments): Nx = int(Nxs[i]) Ny = int(Nys[i]) Nx = Nx+1 Ny = Ny+1 if i < 4: print('%dx%d' % (Nx-1, Ny-1)) print(' to the manufactured solution:') print(" u - um = min/max (%e, %e)" % (numpy.min(u[i] - um[i]), numpy.max(u[i] - um[i]))) print(" v - vm = min/max (%e, %e)" % (numpy.min(v[i] - vm[i]), numpy.max(v[i] - vm[i]))) print(' to daetools:') print(" u - u_dt = min/max (%e, %e)" % (numpy.min(u[i] - u_dt[i]), numpy.max(u[i] - u_dt[i]))) print(" v - v_dt = min/max (%e, %e)" % (numpy.min(v[i] - v_dt[i]), numpy.max(v[i] - v_dt[i]))) print(" um - um_dt = min/max (%e, %e)" % (numpy.min(um[i] - um_dt[i]), numpy.max(um[i] - um_dt[i]))) print(" vm - vm_dt = min/max (%e, %e)" % (numpy.min(vm[i] - vm_dt[i]), numpy.max(vm[i] - vm_dt[i]))) Eu[i] = numpy.sqrt((1.0/(Nx*Ny)) * numpy.sum((u[i]-u_dt[i])**2)) Ev[i] = numpy.sqrt((1.0/(Nx*Ny)) * numpy.sum((v[i]-v_dt[i])**2)) Eum[i] = numpy.sqrt((1.0/(Nx*Ny)) * numpy.sum((um[i]-um_dt[i])**2)) Evm[i] = numpy.sqrt((1.0/(Nx*Ny)) * numpy.sum((vm[i]-vm_dt[i])**2)) print(' Eu(u-u_dt) =', Eu) print(' Ev(v-v_dt) =', Ev) print(' Eum(um-um_dt) =', Eum) print(' Evm(vm-vm_dt) =', Evm) def run(**kwargs): Nxs = numpy.array(Nxs_) Nys = numpy.array(Nys_) n = len(Nxs) Lx = 0.8 hs = Lx / Nxs # It's similar in y direction Eu = numpy.zeros(n) Ev = numpy.zeros(n) Cu = numpy.zeros(n) Cv = numpy.zeros(n) pu = numpy.zeros(n) pv = numpy.zeros(n) E2 = numpy.zeros(n) u_dt = [0]*n v_dt = [0]*n um_dt = [0]*n vm_dt = [0]*n # The normalised global errors if CompareToDAEToolsAndPlotDAETools: for i in range(n): Nx = int(Nxs[i]) Ny = int(Nys[i]) times, u, v, um, vm, simulation = simulate(Nx, Ny, **kwargs) Eu[i] = numpy.sqrt((1.0/(Nx*Ny)) * numpy.sum((u-um)**2)) Ev[i] = numpy.sqrt((1.0/(Nx*Ny)) * numpy.sum((v-vm)**2)) u_dt[i] = u v_dt[i] = v um_dt[i] = um vm_dt[i] = vm u_cs, v_cs, um_cs, vm_cs, Eu_cs, Ev_cs = loadOpenCSData() compareResults(u_cs, v_cs, um_cs, vm_cs, u_dt, v_dt, um_dt, vm_dt) else: u, v, um, vm, Eu, Ev = loadOpenCSData() # Order of accuracy for i,Nx in enumerate(Nxs): if i == 0: pu[i] = 0 pv[i] = 0 Cu[i] = 0 Cv[i] = 0 else: pu[i] = numpy.log(Eu[i]/Eu[i-1]) / numpy.log(hs[i]/hs[i-1]) pv[i] = numpy.log(Ev[i]/Ev[i-1]) / numpy.log(hs[i]/hs[i-1]) Cu[i] = Eu[i] / hs[i]**pu[i] Cv[i] = Ev[i] / hs[i]**pv[i] C2u = 0.030 # constant for the second order slope line (to get close to the actual line) C2v = 0.075 # constant for the second order slope line (to get close to the actual line) Eu2 = C2u * hs**2 # Eu for the second order slope Ev2 = C2v * hs**2 # Ev for the second order slope print('----------------') print('Eu =', Eu) print('Ev =', Ev) print('Eu2 =', Eu2) print('Ev2 =', Ev2) print('pu =', pu) print('pv =', pv) fontsize = 14 fontsize_legend = 11 fig = plt.figure(figsize=(10,8), facecolor='white') grids = ['%dx%d' % (Nx,Ny) for (Nx,Ny) in zip(Nxs,Nys)] grids = ','.join(grids) #fig.canvas.setWindowTitle('The Normalised global errors and the Orders of accuracy (grids: %s) (cv_4)' % grids) ax = plt.subplot(221) plt.figure(1, facecolor='white') plt.loglog(hs, Eu, 'ro', label='Eu(h)') plt.loglog(hs, Eu2, 'b-', label='2nd order slope') plt.xlabel('h', fontsize=fontsize) plt.ylabel('||Eu||', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) #plt.xlim((0.0, 0.09)) ax = plt.subplot(222) plt.figure(1, facecolor='white') plt.semilogx(hs[1:], pu[1:], 'rs-', label='Order of Accuracy (pu)') plt.xlabel('h', fontsize=fontsize) plt.ylabel('pu', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) #plt.xlim((0.0, 0.09)) #plt.ylim((1.94, 2.02)) ax = plt.subplot(223) plt.figure(1, facecolor='white') plt.loglog(hs, Ev, 'ro', label='Ev(h)') plt.loglog(hs, Ev2, 'b-', label='2nd order slope') plt.xlabel('h', fontsize=fontsize) plt.ylabel('||Ev||', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) #plt.xlim((0.0, 0.09)) ax = plt.subplot(224) plt.figure(1, facecolor='white') plt.semilogx(hs[1:], pv[1:], 'rs-', label='Order of Accuracy (pv)') plt.xlabel('h', fontsize=fontsize) plt.ylabel('pv', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) #plt.xlim((0.0, 0.09)) #plt.ylim((1.94, 2.02)) plt.tight_layout() plt.show() if __name__ == "__main__": if len(sys.argv) == 1: CompareToDAEToolsAndPlotDAETools = False elif len(sys.argv) == 2: cmp = sys.argv[1] if cmp == 'opencs': # Loads OpenCS data and plott OpenCS plot. CompareToDAEToolsAndPlotDAETools = False else: # Runs DAE tools model, compares to OpenCS data and plots DAE Tools plot. CompareToDAEToolsAndPlotDAETools = True else: print('Usage: python dae_example_4_cv_plots.py [opencs|comparison]') sys.exit() run()