#! /usr/bin/env python # -*- coding: utf-8 -*- """ Python implementation of an interative PDE solver. The original work was an Octave (Matlab) code, at http://perso.crans.org/besson/publis/PDE_09_2014/pde_09_2014.m. @date: Mon Feb 05 16:15:21 2015. @author: Lilian Besson for CS101 course at Mahindra Ecole Centrale 2015. @licence: GNU Public License version 3. """ from pylab import * clf() # If we want to print symetrically from -1 to 1, else 0 to 1. printSym = False # If we want to use local Frechet average windows. printFrechet = True # Defining the constants k = 1. x0, x1 = 0.0, 1.0 lx = x1 - x0 n = 500 # FIXME change that value for n h = float(lx) / n dt = 5e-4 c = h * dt * sqrt(k) # jmax is the number of executing steps of the numerical simulation jmax = 10000 # FIXME change that value for jmax jskip = 500 # Smoothing the initial profile: useless ? beta = 1.0 #beta = 0.05 # FIXME? # Initializing arrays if printSym: x = arange(-x1, x1, h/2.0) else: x = arange(0.0, x1, h) xMinus = arange(0.0, x1, h) print "\nSize of x is:", len(x) # First initial profile: affine # This is the regular answer proved by the method of caracterists # It will move with a positive translation on the y axis #y0 = beta * (1.0 - x) # FIXME do not use that one # Second initial profile: parabolic # A small pick will grow from the middle point 0 y0 = beta * (1.0 - (x ** 2)) N = len(y0) # Two tools, for the left and right scheme (one scalar version, one componentwise version) iMinusOne = lambda j: max(0, j-1) iPlusOne = lambda j: min(N-1, j+1) iMinusOneV = lambda v: [max(0, j-1) for j in v] iPlusOneV = lambda v: [min(N-1, j+1) for j in v] # Width of the local Frechet average window frechetWindow = 10 # FIXME put 100 again after hFrechet = h*frechetWindow # Size of the local Frechet average nFrechet = round(float((N-1)) / frechetWindow) # Initialisation of the matrix used to compute local Frechet average print "Shape of the Frechet matrix will be", N*N valueMatFrechet = 1.0 / frechetWindow # Thanks to a sparse matrix (we win a LOT of memory) # (Reference is http://docs.scipy.org/doc/scipy/reference/sparse.html) matFrechet = matrix(zeros((N/frechetWindow, N))) # FIXME do something more memory and time efficient to initialize that sparse matrix for i in arange(0, N): for j in filter(lambda k: k == (i / frechetWindow), arange(0, N/frechetWindow)): # print " For i = ", i, " and j =", j # print " Non zero value on matFrechet!" matFrechet[j, i] = valueMatFrechet print "Shape of matFrechet:", size(matFrechet) # The Octave command is sparse(jFrechet, iFrechet, 1/frechet) # jColumn, iColumn, value from scipy import sparse mSparse = sparse.dok_matrix(matFrechet) del matFrechet # We delete the non-sparse matrix: it is useless now print mSparse print mSparse.shape if printSym: xFrechet = arange(-x1, x1, hFrechet) else: xFrechet = arange(x0, x1, hFrechet) print "\nSize of x is:", len(x) # %% Start the simulation! matricialOperations = True # FIXME the componentwise operations seems to not work hold(False) plt.ion() # FIXME (from http://stackoverflow.com/questions/4098131/update-a-plot-in-matplotlib) nameOfSchemes = ['left', 'center', 'right'] # Change here the kind of numerical scheme to be used, between 0, 1, or 2 for schemaN in [0]: schema = nameOfSchemes[schemaN] # Initialisation for the simulation for that scheme clf() fig = figure(schemaN) ax = fig.add_subplot(111) # Init matrix y (column vector) y = matrix(y0).T print 'Shape of the vector y:', shape(y) # Init yFrechet yFrechet = mSparse * y print 'Shape of the second vector yFrechet (for the local Frechet average) :', shape(yFrechet) # Option for the plotting if printFrechet: # We print the local Frechet averages line1, = ax.plot(xFrechet, yFrechet, '-r') # line1.set_ydata(y) ylabel('Frechet average position y of the metal surface, \nyFrechet > 0 (adimensional), with {} values.'.format(nFrechet)) if printSym: xlabel('Position x, -1 < xFrechet < 1 (adimensional), with {} points.'.format(nFrechet)) else: xlabel('Position x, 0 < xFrechet < 1 (adimensional), with {} points.'.format(nFrechet)) else: # We print the entire vector line1, = ax.plot(x, y, '-r') ylabel('Position y of the metal surface, \ny > 0 (adimensional), with {} valeurs.'.format(N)) if printSym: xlabel('Position x, -1 < x < 1 (adimensional), with {} points.'.format(n)) else: xlabel('Position x, 0 < x < 1 (adimensional), with {} points.'.format(n)) grid('on') fig.canvas.draw() # # FIXME if you want to be able to pause in the beginning of the simulation # print raw_input("Can I start the simulation? (with {} steps) ?".format(jmax)) # pause(10) # Starting of the simulation, with jmax steps: for j in xrange(0, jmax + 1): # print "Step j={} for the simulation.".format(j) # Too time consuming if matricialOperations: # We update that vector y (with one matricial computation) if schema == 'left': # 1) Left scheme y += c / sqrt(h**2 + np.power(y - y[iMinusOneV(arange(0, N))], 2)) elif schema == 'center': # 2) Center scheme y += c / sqrt(h**2 + np.power(y[iPlusOneV(arange(0, N))] - y, 2)) elif schema == 'right': # 3) Right scheme y += c / sqrt(h**2 + 0.25 * np.power(y[iPlusOneV(arange(0, N))] - y[iMinusOneV(arange(0, N))], 2)) else: pass else: # Updating y with slow for loops instead of vector operations if schema == 'left': # 1) Left scheme for i in range(0, N): y[i] = y[i] + c / (sqrt(h**2 + (y[i] - y[iMinusOne(i)])**2)) elif schema == 'right': # 2) Right scheme for i in range(0, N): y[i] = y[i] + c / (sqrt(h**2 + (y[iPlusOne(i)] - y[i])**2)) elif schema == 'center': # 3) Center scheme for i in range(0, N): y[i] = y[i] + c / (sqrt(h**2 + 0.25*( y[iPlusOne(i)] - y[iMinusOne(i)] )**2)) else: pass # End of updates for y # Plotting the y vector (from time to time) if (j % jskip) == 0: tt = 'Simulation of the evolution of the metal surface\n in respect to the time (step j={}), with a {} schema \n(C) Generated by Python, code by Lilian Besson and equations by Nicolas Vernier, 2015.\n'.format(j, schema) print j title(tt) if printFrechet: # Smoothing of the approximative solution y (based on a NON-PROVEN local Frechet average) yFrechet = mSparse * y # line1, = ax.plot(xFrechet, yFrechet, '-r') ylim(0.99*min(yFrechet), 1.01*max(yFrechet)) line1.set_ydata(yFrechet) else: # line1, = ax.plot(x, y, '+g') line1.set_ydata(y) fig.show() fig.canvas.draw() # We save the figure: if (j % 10*jskip) == 0: fig.savefig('plot_{}__{}.png'.format(schema, j)) # show() # pause(10) # FIXME print '\n\nEnd of the simulation for the {} scheme.\n'.format(schema) # http://www.mathworks.in/help/matlab/ref/saveas.html fig.savefig('plot_{}__last.png'.format(schema)) fig.savefig('plot_{}__last.svg'.format(schema)) fig.savefig('plot_{}__last.pdf'.format(schema)) pause(5) print '\n\nEnd of the numerical simulation.\nProgram written in Python 2.7+, with NumPy and MatPlotLib, by Lilian Besson (ENS Cachan).\nWritten for a collaborative work with Nicolas Vernier (Université Paris-Sud), at Mahindra École Centrale, (C) 2015.\n' # pause();