#! 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 02 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 = 12400 # FIXME change that value for n h = lx / n dt = 5e-5 c = h * dt * sqrt(k) jmax = 10000 # FIXME change that value for jmax (the number of executing steps of the numerical simulation) jskip = 500 # Smoothing the initial profile 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 translation #y0 = beta * (1.0 - x) # FIXME do not use that one # First initial profile: parabolic # A small pick will grow from the middle point 0 y0 = beta * (1.0 - (x ** 2)) N = len(y0) # Two tools, FIXME useless? iMinusOne = lambda j: max(0, j-1) iPlusOne = lambda j: min(N-1, j+1) #iMinusOneV = lambda v: map(lambda j: max(1, j-1), v) iMinusOneV = lambda v: [ max(0, j-1) for j in v ] #iPlusOneV = lambda v: map(lambda j: min(N, j+1), v) iPlusOneV = lambda v: [ min(N-1, j+1) for j in v ] # Width of the local Frechet average window frechetWindow = 200 # 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) # (Cf. http://docs.scipy.org/doc/scipy/reference/sparse.html) matFrechet = matrix(zeros((N/frechetWindow, N))) # FIXME do something more memory and time efficient for i in arange(0, N): # print "Start with i =", i # for j in arange(0, N/frechetWindow): for j in filter(lambda k: k == (i / frechetWindow), arange(0, N/frechetWindow)): # if j == (i / frechetWindow): # print " For i = ", i, " and j =", j matFrechet[j, i] = valueMatFrechet # print " Non zero value on matFrechet!" # FIXME rewrite correctly! We can do this better! print "Shape of matFrechet:", size(matFrechet) from scipy import sparse mSparse = sparse.dok_matrix(matFrechet) del matFrechet # FIXME ? print mSparse print mSparse.shape # The Octave command is sparse(jFrechet, iFrechet, 1/frechet) # jColumn, iColumn, value # %% if printSym: # xFrechet = -x1 : hFrechet : (x1 - hFrechet) xFrechet = arange(-x1, x1, hFrechet) else: # xFrechet = x0 : hFrechet : (x1 - hFrechet) xFrechet = arange(x0, x1, hFrechet) print "\nSize of x is:", len(x) # %% Start the simulation! matricialOperations = True # FIXME # pause() hold(False) # hold(True) 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 de la simulation pour ce schéma clf() fig = figure(schemaN) # Init y y = matrix(y0).T # FIXME matrix? 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 HandlePlot = plot(xFrechet, yFrechet, '-r') #, 'YDataSource', 'yFrechet') # HandlePlot[0].set_ydata(y) fig.ylabel('Frechet average position y of the metal surface, yFrechet > 0 (adimensional), with {} values.'.format(nFrechet)) if printSym: fig.xlabel('Position x, -1 < xFrechet < 1 (adimensional), with {} points.'.format(nFrechet)) else: fig.xlabel('Position x, 0 < xFrechet < 1 (adimensional), with {} points.'.format(nFrechet)) else: # We print the entire vector HandlePlot = plot(x, y, '-r') #, 'YDataSource', 'y') fig.ylabel('Position y of the metal surface, y > 0 (adimensional), with {} valeurs.'.format(N)) if printSym: fig.xlabel('Position x, -1 < x < 1 (adimensional), with {} points.'.format(n)) else: fig.xlabel('Position x, 0 < x < 1 (adimensional), with {} points.'.format(n)) grid('on') # 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) 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) Schéma gauche for i in range(0, N): y[i] = y[i] + c / (sqrt(h**2 + ( y[i] - y[iMinusOne(i)] )**2)) elif schema == 'right': # 2) Schéma droit for i in range(0, N): y[i] = y[i] + c / (sqrt(h**2 + ( y[iPlusOne(i)] - y[i] )**2)) elif schema == 'center': # 3) Schéma centré 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 # Plotting the y vector (from time to time) if (j % jskip) == 0: tt = 'Simulation of the evolution of the metal surface in respect to the time (step j={}), with a schema {} \n(C) Generate by Python, code by Lilian Besson and equations by Nicolas Vernier, 2015.\n'.format(j, schema) print j fig.title(tt) if printFrechet: # Smoothing of the approximative solution y (based on a NON-PROVEN local Frechet average) yFrechet = mSparse * y # set(HandlePlot, 'YData', yFrechet) HandlePlot = plot(xFrechet, yFrechet, '-r') else: # set(HandlePlot, 'YData', y) HandlePlot = plot(x, y, '+g') show() # We save the figure: if (j % 10*jskip) == 0: fig.savefig('plot_{}__{}.png'.format(schema, j)) fig.show() # fig.draw() # 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();