From 248c5130de9e8451d49480a6ba94c5511fd14e87 Mon Sep 17 00:00:00 2001 From: Thibaut Lunet <thibaut.lunet@unige.ch> Date: Mon, 19 Mar 2018 00:34:01 +0100 Subject: [PATCH] TL : added some first example scripts --- examples/numpy/csv-read-interpol-plot.py | 50 ++++ examples/numpy/data.csv | 302 +++++++++++++++++++++++ examples/scipy/jacobMultiDim.py | 233 +++++++++++++++++ 3 files changed, 585 insertions(+) create mode 100644 examples/numpy/csv-read-interpol-plot.py create mode 100644 examples/numpy/data.csv create mode 100644 examples/scipy/jacobMultiDim.py diff --git a/examples/numpy/csv-read-interpol-plot.py b/examples/numpy/csv-read-interpol-plot.py new file mode 100644 index 0000000..d887dce --- /dev/null +++ b/examples/numpy/csv-read-interpol-plot.py @@ -0,0 +1,50 @@ +# -*- 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 diff --git a/examples/numpy/data.csv b/examples/numpy/data.csv new file mode 100644 index 0000000..d6bb50e --- /dev/null +++ b/examples/numpy/data.csv @@ -0,0 +1,302 @@ +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 diff --git a/examples/scipy/jacobMultiDim.py b/examples/scipy/jacobMultiDim.py new file mode 100644 index 0000000..c3fcc80 --- /dev/null +++ b/examples/scipy/jacobMultiDim.py @@ -0,0 +1,233 @@ +#!/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') -- GitLab