#!/usr/bin/env python # -*- coding: utf-8 -*- """ *********************************************************************************** tutorial_cv_2.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. References: 1. G. Tryggvason. Method of Manufactured Solutions, Lecture 33: Predictivity-I, 2011. `PDF <http://www3.nd.edu/~gtryggva/CFD-Course/2011-Lecture-33.pdf>`_ 2. K. Salari and P. Knupp. Code Verification by the Method of Manufactured Solutions. SAND2000 – 1444 (2000). `doi:10.2172/759450 <https://doi.org/10.2172/759450>`_ 3. P.J. Roache. Fundamentals of Verification and Validation. Hermosa, 2009. `ISBN-10:0913478121 <http://www.isbnsearch.org/isbn/0913478121>`_ The procedure for the *transient convection-diffusion* (Burger's) equation: .. code-block:: none L(f) = df/dt + f*df/dx - D*d2f/dx2 = 0 is the following: 1. Pick a function (q, the manufactured solution): .. code-block:: none q = A + sin(x + Ct) 2. Compute the new source term (g) for the original problem: .. code-block:: none g = dq/dt + q*dq/dx - D*d2q/dx2 3. Solve the original problem with the new source term: .. code-block:: none df/dt + f*df/dx - D*d2f/dx2 = g Since L(f) = g and g = L(q), consequently we have: f = q. Therefore, the computed numerical solution (f) should be equal to the manufactured one (q). The terms in the source g term are: .. code-block:: none dq/dt = C * cos(x + C*t) dq/dx = cos(x + C*t) d2q/dx2 = -sin(x + C*t) The equations are solved for Dirichlet boundary conditions: .. code-block:: none f(x=0) = q(x=0) = A + sin(0 + Ct) f(x=2pi) = q(x=2pi) = A + sin(2pi + Ct) Numerical vs. manufactured solution plot (no. elements = 60, t = 1.0s): .. image:: _static/tutorial_cv_2-results.png :width: 500px The normalised global errors and the order of accuracy plots (no. elements = [60, 90, 120, 150], t = 1.0s): .. image:: _static/tutorial_cv_2-results2.png :width: 800px """ import sys, numpy from time import localtime, strftime from daetools.pyDAE import * import matplotlib.pyplot as plt 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.A = 1.0 self.C = 1.0 self.D = 0.05 self.f = daeVariable("f", no_t, self, "", [self.x]) self.q = daeVariable("q", no_t, self, "", [self.x]) def DeclareEquations(self): daeModel.DeclareEquations(self) # Create some auxiliary functions to make equations more readable A = self.A C = self.C D = self.D t = Time() f = lambda x: self.f(x) df_dt = lambda x: dt(self.f(x)) df_dx = lambda x: d(self.f(x), self.x, eCFDM) d2f_dx2 = lambda x: d2(self.f(x), self.x, eCFDM) q = lambda x: A + numpy.sin(x() + C*t) dq_dt = lambda x: C * numpy.cos(x() + C*t) dq_dx = lambda x: numpy.cos(x() + C*t) d2q_dx2 = lambda x: -numpy.sin(x() + C*t) g = lambda x: dq_dt(x) + q(x) * dq_dx(x) - D * d2q_dx2(x) # Numerical solution eq = self.CreateEquation("f", "Numerical solution") x = eq.DistributeOnDomain(self.x, eOpenOpen) eq.Residual = df_dt(x) + f(x) * df_dx(x) - D * d2f_dx2(x) - g(x) eq.CheckUnitsConsistency = False eq = self.CreateEquation("f(0)", "Numerical solution") x = eq.DistributeOnDomain(self.x, eLowerBound) eq.Residual = f(x) - q(x) eq.CheckUnitsConsistency = False eq = self.CreateEquation("f(2pi)", "Numerical solution") x = eq.DistributeOnDomain(self.x, eUpperBound) eq.Residual = f(x) - q(x) eq.CheckUnitsConsistency = False # Manufactured solution eq = self.CreateEquation("q", "Manufactured solution") x = eq.DistributeOnDomain(self.x, eClosedClosed) eq.Residual = self.q(x) - q(x) eq.CheckUnitsConsistency = False class simTutorial(daeSimulation): def __init__(self, Nx): daeSimulation.__init__(self) self.m = modTutorial("tutorial_cv_2(%d)" % Nx) self.m.Description = __doc__ self.Nx = Nx def SetUpParametersAndDomains(self): self.m.x.CreateStructuredGrid(self.Nx, 0, 2*numpy.pi) def SetUpVariables(self): Nx = self.m.x.NumberOfPoints xp = self.m.x.Points for x in range(1, Nx-1): self.m.f.SetInitialCondition(x, self.m.A + numpy.sin(xp[x])) # 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 # Enable reporting of all variables simulation.m.SetReportingOn(True) # Enable reporting of time derivatives for all reported variables simulation.ReportTimeDerivatives = True # Set the time horizon and the reporting interval simulation.ReportingInterval = 0.05 simulation.TimeHorizon = 1 # 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 fvar = results[simulation.m.Name + '.f'] qvar = results[simulation.m.Name + '.q'] times = fvar.TimeValues q = qvar.Values[-1, :] # 2D array [t,x] f = fvar.Values[-1, :] # 2D array [t,x] #print(times,f,q) return times,f,q def run(**kwargs): Nxs = numpy.array([60, 90, 120, 150]) n = len(Nxs) L = 2*numpy.pi hs = L / Nxs E = numpy.zeros(n) C = numpy.zeros(n) p = numpy.zeros(n) E2 = numpy.zeros(n) # The normalised global errors for i,Nx in enumerate(Nxs): times, numerical_sol, manufactured_sol = simulate(int(Nx)) E[i] = numpy.sqrt((1.0/Nx) * numpy.sum((numerical_sol-manufactured_sol)**2)) # Order of accuracy for i,Nx in enumerate(Nxs): p[i] = numpy.log(E[i]/E[i-1]) / numpy.log(hs[i]/hs[i-1]) C[i] = E[i] / hs[i]**p[i] C2 = 0.09 # constant for the second order slope line (to get close to the actual line) E2 = C2 * hs**2 # E for the second order slope fontsize = 14 fontsize_legend = 11 fig = plt.figure(figsize=(10,4), facecolor='white') fig.canvas.set_window_title('The Normalised global errors and the Orders of accuracy (Nelems = %s) (cv_2)' % Nxs.tolist()) ax = plt.subplot(121) plt.figure(1, facecolor='white') plt.loglog(hs, E, 'ro', label='E(h)') plt.loglog(hs, E2, 'b-', label='2nd order slope') plt.xlabel('h', fontsize=fontsize) plt.ylabel('||E||', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) plt.xlim((0.04, 0.11)) ax = plt.subplot(122) plt.figure(1, facecolor='white') plt.semilogx(hs[1:], p[1:], 'rs-', label='Order of Accuracy (p)') plt.xlabel('h', fontsize=fontsize) plt.ylabel('p', fontsize=fontsize) plt.legend(fontsize=fontsize_legend) plt.xlim((0.04, 0.075)) plt.ylim((2.0, 2.1)) plt.tight_layout() plt.show() if __name__ == "__main__": run()