#!/usr/bin/env python # -*- coding: utf-8 -*- """ *********************************************************************************** tutorial_opencs_dae_6_cv.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__ = """ Code verification using the Method of Manufactured Solutions. It represents re-implementation of DAE Tools code verifications tutorial CV-3 - the 1-D transient convection-diffusion (Burger's) equation. This example uses Neumann boundary conditions. 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>`_ Note: csSimulate cannot be called multiple times from python (multiple MPI init and finalize are not allowed). Therefore, the simulation are performed using the csSimulator binary (it has to be in the PATH for this to work). """ import os, sys, json, itertools, numpy, pandas from daetools.solvers.opencs import pyOpenCS 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 import matplotlib.pyplot as plt class Burgers_1D: def __init__(self, Nx_): self.Nx = Nx_ + 1 Nx = self.Nx self.A = 1.0 self.C = 1.0 self.D = 0.05 self.x0 = 0.0 self.x1 = 2*numpy.pi self.dx = (self.x1-self.x0) / (Nx-1) self.x_domain = [] for x in range(self.Nx): self.x_domain.append(self.x0 + x * self.dx) self.f_start_index = 0*Nx self.q_start_index = 1*Nx self.Nequations = 2*Nx def GetInitialConditions(self): # Use numpy array so that setting f_0 and q_0 changes the original values uv0 = numpy.zeros(self.Nequations) f_0 = uv0[self.f_start_index : self.q_start_index] q_0 = uv0[self.q_start_index : self.Nequations] for ix in range(self.Nx): xp = self.x_domain[ix] f_0[ix] = self.A + numpy.sin(xp) q_0[ix] = self.A + numpy.sin(xp) return uv0.tolist() def GetVariableNames(self): return ['f(%d)'%x for x in range(self.Nx)] + ['q(%d)'%x for x in range(self.Nx)] def CreateEquations(self, y, dydt, currentTime): # 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 f_values = y [self.f_start_index : self.q_start_index] q_values = y [self.q_start_index : self.Nequations] dfdt_values = dydt[self.f_start_index : self.q_start_index] dx = self.dx Nx = self.Nx xp = self.x_domain A = self.A C = self.C D = self.D t = currentTime # Math. functions sin = pyOpenCS.sin cos = pyOpenCS.cos # Variables f = lambda x: f_values[x] qm = lambda x: q_values[x] df_dt = lambda x: dfdt_values[x] def df_dx(x): if(x == 0): # left f0 = f(0) f1 = f(1) f2 = f(2) return (-3*f0 + 4*f1 - f2) / (2*dx) elif(x == Nx-1): # right fn = f(Nx-1) fn1 = f(Nx-1-1) fn2 = f(Nx-1-2) return (3*fn - 4*fn1 + fn2) / (2*dx) else: f1 = f(x+1) f2 = f(x-1) return (f1 - f2) / (2*dx) def d2f_dx2(x): # Interior points fi1 = f(x+1) fi = f(x) fi2 = f(x-1) return (fi1 - 2*fi + fi2) / (dx*dx) q = lambda x: A + sin(xp[x] + C*t) dq_dt = lambda x: C * cos(xp[x] + C*t) dq_dx = lambda x: cos(xp[x] + C*t) d2q_dx2 = lambda x: -sin(xp[x] + C*t) g = lambda x: dq_dt(x) + q(x) * dq_dx(x) - D * d2q_dx2(x) eq = 0 equations = [None] * self.Nequations # Component u: for x in range(Nx): if(x == 0): # BCs: Neumann equations[eq] = df_dx(x) - dq_dx(x) elif(x == Nx-1): # BCs: Neumann equations[eq] = df_dx(x) - dq_dx(x) else: # Interior points equations[eq] = df_dt(x) + f(x) * df_dx(x) - D * d2f_dx2(x) - g(x) eq += 1 # Component v: for x in range(Nx): equations[eq] = qm(x) - q(x) eq += 1 return equations def loadOpenCSData(Nx, resultsDirectory = 'tutorial_opencs_dae_6_cv'): csvFilename = os.path.join(resultsDirectory, '%d' % Nx, '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) Nx = Nx+1 print('Loading %d data...' % (Nx-1)) f_start = 1 # The first value is Time (skip it) f_end = f_start + Nx q_start = f_end q_end = q_start + Nx data = df.values[-1] # take the last time point f = data[f_start:f_end] qm = data[q_start:q_end] return f, qm def simulate(Nx): # Instantiate the model being simulated. model = Burgers_1D(Nx) # 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-6) # 2. Specify the OpenCS model. mb.ModelEquations = model.CreateEquations(mb.Variables, mb.TimeDerivatives, mb.Time) # 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'] = '%d' % Nx options['Simulation']['TimeHorizon'] = 1.0 options['Simulation']['ReportingInterval'] = 0.05 options['Solver']['Parameters']['RelativeTolerance'] = 1e-5 options['LinearSolver'] = { "Library": "Trilinos", "Name": "Amesos_Klu", "PrintInfo": False, "Parameters": {} } 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 inputFilesDirectory = 'tutorial_opencs_dae_6_cv' cs_models = mb.PartitionSystem(Npe, graphPartitioner) csModelBuilder_t.ExportModels(cs_models, inputFilesDirectory, mb.SimulationOptions) print("Single CPU OpenCS model generated successfully!") # 4. Run simulation using the exported model from the specified directory. # Cannot call csSimulate multiple times from python (multiple MPI init and finalize are not allowed). # Therefore, simulate it using the csSimulator binary (it has to be in the PATH for this to work). #csSimulate(inputFilesDirectory) os.system('csSimulator %s' % inputFilesDirectory) f, q = loadOpenCSData(Nx) return 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): 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.18 # 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.setWindowTitle('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.03)) plt.tight_layout() plt.show() if __name__ == "__main__": run()