TL : added some first example scripts

examples/numpy/csv-read-interpol-plot.py

+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jun 25 17:28:03 2015
+
+Python script for post-processing data stored in a csv file
+
+@author: lunet
+"""
+# Importing all required libraries
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Get data from the csv file
+data = np.genfromtxt('data.csv',dtype=float,delimiter=',',names=True)
+
+# Getting the fields
+times = data['Time']
+pressures = data['Pressure']
+velocities = data['Velocity']
+
+# Interpoling data to suppress measurement error
+polyDeg = 25
+# getting the coefficient of the best polynomial fit
+pressuresPolyCoeff = np.polyfit(times,pressures,polyDeg)
+velocitiesPolyCoeff = np.polyfit(times,velocities,polyDeg)
+
+# Constructing the polynomial with the coefficients
+pressuresPoly=np.poly1d(pressuresPolyCoeff)
+velocitiesPoly=np.poly1d(velocitiesPolyCoeff)
+
+# Evaluate the polynomial at each time
+pressuresInter = pressuresPoly(times)
+velocitiesInter = velocitiesPoly(times)
+
+# Plotting all data
+plt.figure()
+plt.subplot(211)
+plt.plot(times,pressures,'b',label="Measurements")
+plt.plot(times,pressuresInter,'r',label="Interpolation")
+plt.xlabel('Time')
+plt.ylabel('Pressure')
+plt.title('Comparison between measured and interpolated fields')
+plt.legend(loc=4)
+plt.subplot(212)
+plt.plot(times,velocities,'b',label="Measurements")
+plt.plot(times,velocitiesInter,'r',label="Interpolation")
+plt.xlabel('Time')
+plt.ylabel('Velocity')
+plt.legend(loc=1)
+plt.show()
\ No newline at end of file

examples/numpy/data.csv

+Time,Pressure,Velocity
+0,0.3301958987,3.3338040438
+0.1,0.4168124765,3.4942511735
+0.2,0.4199741904,3.3318186404
+0.3,0.3999769621,3.4085712565
+0.4,0.3169812118,3.4545081834
+0.5,0.4054077216,3.3519612524
+0.6,0.3371851687,3.3642938982
+0.7,0.382838647,3.4543539478
+0.8,0.3801052262,3.3437793355
+0.9,0.3977826942,3.4440860511
+1,0.443507798,3.4928140727
+1.1,0.4720666843,3.4535783909
+1.2,0.3358170144,3.4034250161
+1.3,0.456112626,3.3339432065
+1.4,0.3049156874,3.489547021
+1.5,0.3253959296,3.444929609
+1.6,0.4975764383,3.2986570456
+1.7,0.4490807824,3.462904467
+1.8,0.4550250256,3.2995482411
+1.9,0.4621653188,3.3183618565
+2,0.4150395444,3.4281186669
+2.1,0.4658928536,3.4409782783
+2.2,0.315028259,3.3483666398
+2.3,0.3331350024,3.3546070275
+2.4,0.3233213783,3.4520355336
+2.5,0.3558313471,3.3282919806
+2.6,0.3416307801,3.3301428461
+2.7,0.3723347955,3.4426342946
+2.8,0.5007305806,3.4602381571
+2.9,0.3297186657,3.3528532623
+3,0.4595264414,3.3006331728
+3.1,0.3291531727,3.3566324756
+3.2,0.3572903992,3.3810220598
+3.3,0.406259343,3.4303815242
+3.4,0.4810322711,3.4587080233
+3.5,0.4575413236,3.3451688053
+3.6,0.4969042561,3.3905372445
+3.7,0.3194688803,3.3631406513
+3.8,0.4865212982,3.3706099321
+3.9,0.3541417437,3.3646170199
+4,0.4790020519,3.4528077654
+4.1,0.4097893374,3.3994406498
+4.2,0.5190777686,3.4220160089
+4.3,0.3520860929,3.3299352512
+4.4,0.4579296429,3.2958890771
+4.5,0.4927083143,3.336221616
+4.6,0.3521294512,3.4585504677
+4.7,0.4849063093,3.3532872958
+4.8,0.4017033728,3.4284135169
+4.9,0.3786548571,3.3313316614
+5,0.4579963272,3.3832224289
+5.1,0.4595545255,3.4133673863
+5.2,0.457888411,3.3681513746
+5.3,0.4357921244,3.4611317584
+5.4,0.3671121605,3.4155559521
+5.5,0.5124545759,3.4435188654
+5.6,0.3582031451,3.2696983827
+5.7,0.4417886804,3.285811134
+5.8,0.418331507,3.4573031316
+5.9,0.3878851461,3.4077309123
+6,0.4534186309,3.2972311987
+6.1,0.4849541065,3.3683661687
+6.2,0.3839811902,3.4076273183
+6.3,0.3585842407,3.383589391
+6.4,0.3824551819,3.3400535267
+6.5,0.4322180106,3.3351955876
+6.6,0.4075950103,3.4132416921
+6.7,0.5174642781,3.4376871922
+6.8,0.5495609858,3.4103095359
+6.9,0.3659575787,3.3573699466
+7,0.4221904702,3.3107171569
+7.1,0.4941202544,3.3187270817
+7.2,0.4606456562,3.4284094725
+7.3,0.4684463747,3.2859828593
+7.4,0.4076954533,3.2608252809
+7.5,0.4664817227,3.3232559719
+7.6,0.5200890942,3.3499717249
+7.7,0.5251115496,3.3234949684
+7.8,0.5556174632,3.3842764814
+7.9,0.4121317549,3.2350533476
+8,0.5411558718,3.2379630537
+8.1,0.5287391543,3.3374298915
+8.2,0.5611765999,3.3676337712
+8.3,0.5567703226,3.2872965803
+8.4,0.4477245588,3.2991315247
+8.5,0.4147148519,3.2908893173
+8.6,0.4355303299,3.2776401775
+8.7,0.4277769755,3.2357731933
+8.8,0.5646313399,3.2618131671
+8.9,0.4372415258,3.3778659271
+9,0.5471318275,3.2806715536
+9.1,0.4324147216,3.3679025954
+9.2,0.4071256051,3.2451734462
+9.3,0.5457502753,3.2778193218
+9.4,0.5345401931,3.2968031975
+9.5,0.5807317846,3.2058502601
+9.6,0.5540955274,3.3521460809
+9.7,0.5553435776,3.269231731
+9.8,0.5903828029,3.2754617447
+9.9,0.5455639432,3.2156104879
+10,0.5957009573,3.1870846801
+10.1,0.5569343582,3.2280311662
+10.2,0.513267197,3.2793892946
+10.3,0.5315054505,3.2297164584
+10.4,0.5512933708,3.2133862253
+10.5,0.5494022022,3.1637788571
+10.6,0.5315946083,3.3106984163
+10.7,0.5419697761,3.2024889391
+10.8,0.6202755755,3.1918106838
+10.9,0.5365915684,3.2983374316
+11,0.5735807284,3.3243900849
+11.1,0.5615925685,3.1897222621
+11.2,0.5899966799,3.1319208054
+11.3,0.6835816386,3.2706556642
+11.4,0.5759022411,3.2645663884
+11.5,0.5146287744,3.2270708301
+11.6,0.6551129874,3.1692203097
+11.7,0.6210250221,3.1779980088
+11.8,0.619848244,3.2455433509
+11.9,0.580717813,3.1815626667
+12,0.6294333297,3.175601506
+12.1,0.6530877001,3.2207009133
+12.2,0.7332279929,3.2057456424
+12.3,0.5998685929,3.1733357111
+12.4,0.6198231323,3.1241547452
+12.5,0.7689441369,3.019718749
+12.6,0.7768047092,3.0447362768
+12.7,0.6847395543,3.1209192442
+12.8,0.7415889201,3.1391043652
+12.9,0.7833815158,3.0981268006
+13,0.7787413064,3.0883522516
+13.1,0.8951509063,2.9939702169
+13.2,0.7871184347,2.8842741875
+13.3,0.8975281109,3.0296025059
+13.4,0.7957244461,2.9665465882
+13.5,0.9954188386,2.9250997385
+13.6,0.9955352957,2.81127638
+13.7,0.9879790311,2.7542940756
+13.8,1.0459809396,2.7139586551
+13.9,1.1278747597,2.8292976492
+14,1.1450967873,2.752139278
+14.1,1.1773017906,2.6449437314
+14.2,1.2819072876,2.5204500626
+14.3,1.3846800159,2.5328863229
+14.4,1.2620436669,2.448316339
+14.5,1.5064780133,2.4285760144
+14.6,1.4709124339,2.3136485421
+14.7,1.6233117806,2.152866246
+14.8,1.6126686267,2.057825569
+14.9,1.7221960363,2.0011529831
+15,1.9883204921,1.9476558505
+15.1,2.0521147732,1.7162971315
+15.2,2.0226306914,1.7816078069
+15.3,2.1572578487,1.6540647111
+15.4,2.2569498465,1.4642268561
+15.5,2.3432748245,1.4420361512
+15.6,2.3966306635,1.3532678678
+15.7,2.5457892593,1.2913397949
+15.8,2.5725168304,1.244452223
+15.9,2.5388159546,1.2646271752
+16,2.752607866,1.1729153571
+16.1,2.6355285541,1.041542411
+16.2,2.7631273064,1.0700012342
+16.3,2.8729254089,1.0328740974
+16.4,2.856345997,0.8873737385
+16.5,2.9644585469,0.9121833045
+16.6,2.9201572096,0.9270660978
+16.7,2.8626039484,0.8990080633
+16.8,2.8882057832,0.7502550971
+16.9,3.055864953,0.7897890418
+17,2.9209683923,0.7372343491
+17.1,3.0528734185,0.8214419064
+17.2,3.0620259392,0.6630270028
+17.3,3.1025448965,0.7555867387
+17.4,3.045512489,0.6901407345
+17.5,3.1802519311,0.6804170513
+17.6,3.137626082,0.7517297483
+17.7,3.084271575,0.7694649924
+17.8,3.0757982933,0.7204168631
+17.9,3.1890550403,0.6800974591
+18,3.247843294,0.5845964789
+18.1,3.1222543355,0.6654019505
+18.2,3.0766583471,0.7224440195
+18.3,3.1461555153,0.5269357773
+18.4,3.0872163542,0.6938109104
+18.5,3.184346326,0.5703389962
+18.6,3.158652294,0.6687895388
+18.7,3.1097630489,0.5747569612
+18.8,3.3048795413,0.5335158642
+18.9,3.2189873375,0.5645434357
+19,3.2341890086,0.6555727297
+19.1,3.1465317489,0.6176805787
+19.2,3.3139080729,0.6362504578
+19.3,3.2340090579,0.5249942709
+19.4,3.2719991067,0.5849029669
+19.5,3.1978980666,0.4554938878
+19.6,3.2885345554,0.5042731003
+19.7,3.1930726214,0.5485447885
+19.8,3.3076787835,0.4727096319
+19.9,3.1826486567,0.501238585
+20,3.2146050423,0.4661789
+20.1,3.2663946238,0.4464840521
+20.2,3.2719419178,0.5489247596
+20.3,3.2593143102,0.5811836321
+20.4,3.2813683651,0.5318850211
+20.5,3.3834932616,0.4793161074
+20.6,3.378894372,0.5901062574
+20.7,3.2036550234,0.443745261
+20.8,3.230355234,0.502206709
+20.9,3.2684819126,0.5146087268
+21,3.3798570648,0.5031386604
+21.1,3.3956836405,0.5351312777
+21.2,3.3754833132,0.4535007143
+21.3,3.2623279276,0.5446557172
+21.4,3.3092460304,0.5574030214
+21.5,3.2628528574,0.5091666309
+21.6,3.3717998892,0.484558308
+21.7,3.2415856711,0.5209023832
+21.8,3.2838610693,0.5587477889
+21.9,3.2748486315,0.4954176173
+22,3.2711257816,0.5376572564
+22.1,3.4230038538,0.4441655906
+22.2,3.2344665856,0.4856581876
+22.3,3.2961253192,0.4120874855
+22.4,3.3079764548,0.545818655
+22.5,3.429530703,0.3840614827
+22.6,3.2903143634,0.4683460067
+22.7,3.3510455389,0.5387355846
+22.8,3.2841744362,0.5427496237
+22.9,3.3280108335,0.4698360613
+23,3.2936460703,0.4153486908
+23.1,3.3528926652,0.3988565163
+23.2,3.3837372365,0.3544485043
+23.3,3.2813107338,0.3520348277
+23.4,3.4168456909,0.3522403315
+23.5,3.3555930094,0.5285085731
+23.6,3.4010172136,0.5262153523
+23.7,3.295324915,0.51240149
+23.8,3.3864454895,0.4485649008
+23.9,3.4570886942,0.4631315718
+24,3.3727396774,0.5218367252
+24.1,3.2953087754,0.396749482
+24.2,3.3558506294,0.4296639164
+24.3,3.3047085979,0.4637693304
+24.4,3.2704831328,0.5063152976
+24.5,3.3228087966,0.3899608156
+24.6,3.3699547668,0.3889402233
+24.7,3.3856960316,0.5222049005
+24.8,3.3106942826,0.4940120931
+24.9,3.2841069906,0.3378372216
+25,3.4052832622,0.4434933991
+25.1,3.4207024091,0.4700080973
+25.2,3.2866346036,0.4735196187
+25.3,3.3779772465,0.4513054844
+25.4,3.3412365526,0.5173266685
+25.5,3.2965070516,0.339221913
+25.6,3.3838829747,0.3631846351
+25.7,3.4198034533,0.406769805
+25.8,3.407057509,0.4532342728
+25.9,3.4444875691,0.3891645244
+26,3.433469351,0.4458252365
+26.1,3.3518872499,0.4264275946
+26.2,3.4354769312,0.4675635031
+26.3,3.3792975088,0.4206960593
+26.4,3.3832605649,0.4402689784
+26.5,3.4419269238,0.4040746201
+26.6,3.4557976486,0.440775627
+26.7,3.3662930764,0.5121539483
+26.8,3.3990025964,0.4474377416
+26.9,3.2876446955,0.3406723515
+27,3.3832421683,0.4325302672
+27.1,3.4439804168,0.4605083596
+27.2,3.4868296416,0.3218940493
+27.3,3.4115243434,0.4834023408
+27.4,3.4177223693,0.5024619786
+27.5,3.4390578138,0.3801045742
+27.6,3.4111411391,0.4544919225
+27.7,3.3659585085,0.5071535658
+27.8,3.4505015467,0.4675602237
+27.9,3.3594910719,0.4471351661
+28,3.3100503018,0.4989461888
+28.1,3.3634129859,0.3929810082
+28.2,3.3703970791,0.3375825783
+28.3,3.3245667336,0.3316282637
+28.4,3.4618495777,0.3960343105
+28.5,3.4950988336,0.3994945075
+28.6,3.3733787127,0.3578525181
+28.7,3.3463187324,0.3599008148
+28.8,3.4531648082,0.3476803319
+28.9,3.3823018278,0.4112261687
+29,3.3214369814,0.3263771201
+29.1,3.387284491,0.4494437722
+29.2,3.3468382504,0.3530722701
+29.3,3.3041608407,0.4614934473
+29.4,3.3075409295,0.3980132617
+29.5,3.3099612734,0.4925891614
+29.6,3.3525022323,0.3300722077
+29.7,3.432606178,0.3858062324
+29.8,3.4743102942,0.3528050621
+29.9,3.3909901022,0.3528397049
+30,3.4913687717,0.4875921154

examples/scipy/jacobMultiDim.py

+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+"""
+Created on Tue May 23 13:29:54 2017
+
+@author: t.lunet
+"""
+import numpy as np
+import scipy.linalg as spl
+import matplotlib.pyplot as plt
+
+# Matplotlib : do not fill marker in plot
+plt.rcParams['markers.fillstyle'] = 'none'
+
+# Number of point in each direction
+nX = 12
+nY = 12
+nZ = 12
+
+# Velocity in each direction
+cX = 1.0
+cY = 1.0
+cZ = 1.0
+
+# Discretization scheme in each direction
+schemeX = 'U1'
+schemeY = 'U1'
+schemeZ = 'U1'
+
+# Script actions
+plotMatrix = False
+saveFig = False
+
+# Initial velocity fields
+u01D = np.ones(nX)
+u02D = np.ones((nX, nY)).ravel()
+u03D = np.ones((nX, nY, nZ)).ravel()
+
+# Stencils
+stencilU1 = [-1, 1, 0]
+stencilC2 = [-0.5, 0, 0.5]
+stencilU3 = [1./6, -6./6, 3./6, 2./6, 0]
+
+dicoStencil = {'U1': ['$1^{st}$ order Upwind', stencilU1],
+               'C2': ['$2^{nd}$ order Centered', stencilC2],
+               'U3': ['$3^{rd}$ order Upwind', stencilU3]}
+
+stencilXName = dicoStencil[schemeX][0]
+stencilYName = dicoStencil[schemeY][0]
+stencilZName = dicoStencil[schemeZ][0]
+
+stencilX = dicoStencil[schemeX][1]
+stencilY = dicoStencil[schemeY][1]
+stencilZ = dicoStencil[schemeZ][1]
+
+iSX = (len(stencilX)-1)//2
+iSY = (len(stencilY)-1)//2
+iSZ = (len(stencilZ)-1)//2
+
+
+def rhs1D(u):
+    """
+    Define a 1D Right Hand Side (RHS) for the advection operator
+
+    .. math::
+        f(u) = c_x \\frac{\\partial u}{\\partial x}
+
+    Parameter
+    ---------
+    u : numpy.ndarray
+        The 1D vector, to evaluate the space derivative from,
+        of size :math:`n_x`
+
+    Returns
+    -------
+    uEval : numpy.ndarray
+        The 1D derivative of u
+    """
+    u1D = u.reshape((nX))
+    uEval = np.zeros_like(u1D)
+    for i in range(len(stencilX)):
+        uEval -= stencilX[i]*np.roll(u1D, iSX-i, axis=0)*cX
+
+    return uEval.ravel()
+
+
+def rhs2D(u):
+    """
+    Define a 2D Right Hand Side (RHS) for the advection operator
+
+    .. math::
+        f(u) = c_x \\frac{\\partial u}{\\partial x} +
+        c_y \\frac{\\partial u}{\\partial y}
+
+    Parameter
+    ---------
+    u : numpy.ndarray
+        The 2D vector, to evaluate the space derivative from,
+        of size :math:`n_x.n_y`
+
+    Returns
+    -------
+    uEval : numpy.ndarray
+        The 2D derivative of u
+    """
+    u2D = u.reshape((nX, nY))
+    uEval = np.zeros_like(u2D)
+
+    for i in range(len(stencilX)):
+        uEval -= stencilX[i]*np.roll(u2D, iSX-i, axis=0)*cX
+    for i in range(len(stencilY)):
+        uEval -= stencilY[i]*np.roll(u2D, iSY-i, axis=1)*cY
+
+    return uEval.ravel()
+
+
+def rhs3D(u):
+    """
+    Define a 3D Right Hand Side (RHS) for the advection operator
+
+    .. math::
+        f(u) = c_x \\frac{\\partial u}{\\partial x} +
+        c_y \\frac{\\partial u}{\\partial y} +
+        c_z \\frac{\\partial u}{\\partial z}
+
+    Parameter
+    ---------
+    u : numpy.ndarray
+        The 3D vector, to evaluate the space derivative from,
+        of size :math:`n_x.n_y.n_z`
+
+    Returns
+    -------
+    uEval : numpy.ndarray
+        The 3D derivative of u
+    """
+    u3D = u.reshape((nX, nY, nZ))
+    uEval = np.zeros_like(u3D)
+
+    for i in range(len(stencilX)):
+        uEval -= stencilX[i]*np.roll(u3D, iSX-i, axis=0)*cX
+    for i in range(len(stencilY)):
+        uEval -= stencilY[i]*np.roll(u3D, iSY-i, axis=1)*cY
+    for i in range(len(stencilZ)):
+        uEval -= stencilZ[i]*np.roll(u3D, iSZ-i, axis=2)*cZ
+
+    return uEval.ravel()
+
+
+def computeJacobian(u, rhs):
+    """
+    Compute the Jacobi matrix of a given RHS function, by computing the
+    partial derivative with finite difference for every element of the
+    vectorial basis.
+
+    Parameter
+    ---------
+    u : numpy.ndarray
+        The initial vector, where to evaluate
+    """
+    n = u.size
+    uPerturb = u.copy()
+    jacobian = np.empty((n, n))
+    uEval0 = rhs(u)
+    eps = 1e-8
+    for i in range(n):
+        uPerturb[i] += eps
+        jacobian[:, i] = rhs(uPerturb)
+        jacobian[:, i] -= uEval0
+        jacobian[:, i] /= eps
+        uPerturb[i] = u[i]
+    return jacobian
+
+# Compute the Jacobie matrix of each RHS function
+print('Compute 1D Jacobi matrix')
+A1D = computeJacobian(u01D, rhs1D)
+print('Compute 2D Jacobi matrix')
+A2D = computeJacobian(u02D, rhs2D)
+print('Compute 3D Jacobi matrix')
+A3D = computeJacobian(u03D, rhs3D)
+
+# Eventually plot the Jacobi matrices
+if plotMatrix:
+    plt.figure('1D')
+    plt.imshow(A1D)
+    plt.colorbar()
+
+    plt.figure('2D')
+    plt.imshow(A2D)
+    plt.colorbar()
+
+    plt.figure('3D')
+    plt.imshow(A3D)
+    plt.colorbar()
+
+# Compute the eigenvalues of each Jacobi matrix
+print('Compute eigenvalues of 1D Jacobi matrix')
+lam1D = spl.eigvals(A1D)
+print('Compute eigenvalues of 2D Jacobi matrix')
+lam2D = spl.eigvals(A2D)
+print('Compute eigenvalues of 3D Jacobi matrix')
+lam3D = spl.eigvals(A3D)
+print('Done, ploting ...')
+
+# Plot the eigenvalues
+plt.figure('Eigenvalues')
+plt.axvline(0, c='gray', ls='--', lw=1.0)
+plt.axhline(0, c='gray', ls='--', lw=1.0)
+plt.plot(lam1D.real, lam1D.imag, 's', label='1D', ms=12.0, mfc=None)
+plt.plot(lam2D.real, lam2D.imag, 'o', label='2D', ms=8.0, mfc=None)
+plt.plot(lam3D.real, lam3D.imag, '*', label='3D', ms=6.0, mfc=None)
+plt.axis('equal')
+plt.legend()
+plt.xlabel('$\\mathcal{R}(\\lambda)$')  # Use $...$ for latex output
+plt.ylabel('$\\mathcal{I}(\\lambda)$')
+
+# Define figure title
+if stencilXName == stencilYName and stencilYName == stencilZName:
+    stencilName = stencilXName
+else:
+    stencilName = '{}(x), {}(y), {}(z)'.format(
+        stencilXName, stencilYName, stencilXName)
+
+plt.title('Eigenvalues of '+stencilName+' operator (Advection)\n' +
+          'cX={}, cY={}, cZ={}'.format(cX, cY, cZ))
+
+# Eventually save figure
+if saveFig:
+    if schemeX == schemeY and schemeY == schemeZ:
+        scheme = schemeX
+    else:
+        scheme = schemeX+schemeY+schemeZ
+    plt.savefig('eigenvalues{}.pdf'.format(scheme), bbox_inches='tight')