pencil.calc =========== .. py:module:: pencil.calc .. autoapi-nested-parse:: Math functions and further calculations. Submodules ---------- .. toctree:: :maxdepth: 1 /code/sourcePython/pencil/calc/Gaussian_averages/index /code/sourcePython/pencil/calc/Reynolds/index /code/sourcePython/pencil/calc/accuracy/index /code/sourcePython/pencil/calc/aver2h5/index /code/sourcePython/pencil/calc/draglift/index /code/sourcePython/pencil/calc/example_shocktube/index /code/sourcePython/pencil/calc/fill_gaps_in_grid/index /code/sourcePython/pencil/calc/gravity/index /code/sourcePython/pencil/calc/part_to_grid/index /code/sourcePython/pencil/calc/shocktube/index /code/sourcePython/pencil/calc/streamlines/index /code/sourcePython/pencil/calc/tensors/index Attributes ---------- .. autoapisummary:: pencil.calc.pc_print Classes ------- .. autoapisummary:: pencil.calc.Stream pencil.calc.Draglift pencil.calc.Tensors Functions --------- .. autoapisummary:: pencil.calc.part_to_grid pencil.calc.fill_gaps_in_grid pencil.calc.twonorm_accuracy pencil.calc.order_accuracy pencil.calc.twonorm_array pencil.calc.twonorm pencil.calc.maxerror pencil.calc.twonorm_accuracy2D pencil.calc.twonorm_accuracy1D pencil.calc.draglift pencil.calc.find_extrema pencil.calc.grav_profile pencil.calc.tensors_sph pencil.calc.dot pencil.calc.dot2 pencil.calc.cross pencil.calc.div pencil.calc.curl pencil.calc.curl2 pencil.calc.grad pencil.calc.del2 pencil.calc.del6 pencil.calc.fluid_reynolds pencil.calc.magnetic_reynolds pencil.calc.sod pencil.calc.kernel_smooth pencil.calc.gauss_3Dsmooth Package Contents ---------------- .. py:data:: pc_print .. py:class:: Stream(field, params, xx=(0, 0, 0), time=(0, 1), metric=None, splines=None) Bases: :py:obj:`object` Contains the methods and results for the streamline tracing for a field on a grid. Trace a field starting from xx in any rectilinear coordinate system with constant dx, dy and dz and with a given metric. call signature: tracer(field, params, xx=[0, 0, 0], time=[0, 1], metric=None, splines=None): Keyword arguments: *field*: Vector field which is integrated over with shape [n_vars, nz, ny, nx]. Its elements are the components of the field using unnormed unit-coordinate vectors. *params*: Simulation and tracer parameters. *xx*: Starting point of the field line integration with starting time. *time*: Time array for which the tracer is computed. *metric*: Metric function that takes a point [x, y, z] and an array of shape [3, 3] that has the comkponents g_ij. Use 'None' for Cartesian metric. *splines*: Pre-computed spline coefficients for tricubic interpolation. This can speed up the calculations greatly for repeated streamline tracing on the same data. Should be computed using scipy.ndimage.spline_filter with order=3 for each vector component, i.e.: splines = np.array([spline_filter(field[0], order=3), spline_filter(field[1], order=3), spline_filter(field[2], order=3)]) .. py:attribute:: params .. py:attribute:: xx :value: (0, 0, 0) .. py:attribute:: time :value: (0, 1) .. py:attribute:: section_l .. py:attribute:: total_l .. py:attribute:: iterations .. py:attribute:: section_dh .. py:attribute:: total_h :value: 1 .. py:function:: part_to_grid(xp, yp, zp=False, quantity=False, Nbins=[1024, 1024, 1024], sim=False, extent=False, fill_gaps=False) Bins quantity based on position data xp and yp to 1024^2 bins like a histrogram. This method is not using TSC. :param - xp: array of x and y positions :param yp: array of x and y positions :param - zp: specify if 3D run, set False to have zp == 0 :param - quantity: array of same shape as xp and yp, but with quantity to bin, set it False to count number of occurrences/histrogram2d :param - Nbins: number of histrogram bins for each direction. if 2d only the first two entries in Nbins are used :param - sim: to extract extent from by loading ghostzone-free grid :param - extent: [[xmin, xmax],[ymin, ymax]] or set false and instead give a sim set extent manually e.g. if you want to include ghost zones :param - fill_gaps interpolate empty grid cells: Returns: arr, xgrid, ygrid - arr: 2d array with binned values - x-/ygrid: linspace of used x/y grid - zgrid: if zp != False .. rubric:: Example vpx = part_to_grid_2d(pvar.xp, pvar.yp, pvar.vpx), notice that this will execute pc.get_sim() internally to get the extent .. py:function:: fill_gaps_in_grid(array, key='nan', steps=['interpolate', 'extrapolate'], order=1, DEBUG=False) Interpolates nans by nearest neighbor method to fill gaps in arrays. Beware this method does not invoke bondaries correctly! This method does not work with array beeing a vector field! Hence, use componentwise and stitch together manually! :param - array: array with nans inside :param - key: indicates array entries that will be interpolated :param - method: as specified in scipy.interpolate.griddata, i.e. use 'linear' for 3D data :returns: array without nans inside. .. rubric:: Example [[0, 1, 0], [[0, 1, 0,], [1, nan, 1], -interpol-> [1, 2, 1], [0, 1, 0]] [0, 1, 0]] .. py:function:: twonorm_accuracy(simulations, field='ux', strip=0, var_file='ogvar.dat', direction='x', noerr=True, quiet=True) Assessment of accuracy of simulation: Computes the two-norm error of all available simulation, where the simulation with the maximum amount of grid points is used as the correct/reference solution. E.g., for runs with grid points assessment of accuracy of x-component of velocity along the y-direction, for runs with grid points nxgrid = n, 2n, 4n, 8n, compute || u_n - u_0 || = dy \sum\limits_{n=0}^n (sqrt(u_n(x_n)-u_0(x_8n))) for all runs (except for the 8n case, used as reference). Requires that the runs have matching grids, that is, grid refined by a factor of 2m, and grids adjusted so that the grid point overlap (needs ofset if periodic BC is used). call signature: twonorm_accuracy(simulations) Keyword arguments: *simulations* array of simulation names to be included in the computations *field* variable used in accuracy assessment *strip*: index for strip along coordinate *var_file*: name of varfile to read from each sim *direction*: compute two-norm along 'x' or 'y' direction *noerr*: set to false if you want to return an array of maximum error along strip, in addition to the two-norm Returns array of two-norms where the larges array is used as base .. py:function:: order_accuracy(simulations=[], nstrips=0, twonorm_arr=[], field='ux', var_file='ogvar.dat', direction='x') Compute an estimate of the order of accuracy, using two-norms where the finest grid is used as reference solution u_0. Return array of orders of accuracy, where each order p is for computation along one strip. p = log(||u(\Delta x) - u_0|| / ||u(\Delta x/2) - u_0||)/log(2) Keyword arguments: *simulations* array of simulation names to be included in the computations *nstrips* number of strips to include in twonorm_array twonorm is computed along each strip, and the strips must coencide on the each grid used in such a grid refinement study *twonorm_arr* array of computed two-normes, if these are available if empty, routine will compute these *field* variable used in accuracy assessment *varfile*: name of varfile to read from each sim *direction*: compute two-norm along 'x' or 'y' direction .. py:function:: twonorm_array(simulations, nstrips, field='ux', var_file='ogvar.dat', direction='x') Compute twonorm_accuracy along the selected direction, for a number of 'strips'. See twonorm_accuracy for details. call signature: twonorm_array(simulations,nstrips) Keyword arguments: *simulations* array of simulation names to be included in the computations *nstrips* number of strips to include in twonorm_array twonorm is computed along each strip, and the strips must coencide on the each grid used in such a grid refinement study *field* variable used in accuracy assessment *varfile*: name of varfile to read from each sim *direction*: compute two-norm along 'x' or 'y' direction .. py:function:: twonorm(u1, u2, dx) Compute the two-norm error of two spatial vectors u1 and u2. The distance between grid points used in calculation is dx. .. py:function:: maxerror(u1, u2) Find the larges error between values from u1 and u2 (for the same indices). .. py:function:: twonorm_accuracy2D(simulations, field='ur', var_file='ogvar.dat', noerr=True, quiet=True) Assessment of accuracy of simulation: Computes the two-norm error of all available simulation, where the simulation with the maximum amount of grid points is used as the correct/reference solution. E.g., for runs with grid points assessment of accuracy of x-component of velocity along the y-direction, for runs with grid points nxgrid = n, 2n, 4n, 8n, compute || u_n - u_0 || = dy \sum\limits_{n=0}^n (sqrt(u_n(x_n)-u_0(x_8n))) for all runs (except for the 8n case, used as reference). Requires that the runs have matching grids, that is, grid refined by a factor of 2m, and grids adjusted so that the grid point overlap (needs ofset if periodic BC is used). call signature: twonorm_accuracy(simulations) Keyword arguments: *simulations* array of simulation names to be included in the computations *field* variable used in accuracy assessment *varfile*: name of varfile to read from each sim *noerr*: set to false if you want to return an array of maximum error along strip, in addition to the two-norm Returns array of two-norms where the larges array is used as base .. py:function:: twonorm_accuracy1D(simulations, field='ur', strip=1, direction='r', varfile='ogvar.dat', noerr=True, quiet=True) Assessment of accuracy of simulation: Computes the two-norm error of all available simulation, where the simulation with the maximum amount of grid points is used as the correct/reference solution. E.g., for runs with grid points assessment of accuracy of x-component of velocity along the y-direction, for runs with grid points nxgrid = n, 2n, 4n, 8n, compute || u_n - u_0 || = dy \sum\limits_{n=0}^n (sqrt(u_n(x_n)-u_0(x_8n))) for all runs (except for the 8n case, used as reference). Requires that the runs have matching grids, that is, grid refined by a factor of 2m, and grids adjusted so that the grid point overlap (needs ofset if periodic BC is used). call signature: twonorm_accuracy(simulations) Keyword arguments: *simulations* array of simulation names to be included in the computations *field* variable used in accuracy assessment *varfile*: name of varfile to read from each sim *noerr*: set to false if you want to return an array of maximum error along strip, in addition to the two-norm Returns array of two-norms where the larges array is used as base .. py:function:: draglift(simulations, sortby='dx', **kwargs) Compute mean drag coefficient, rms lift coefficient, and Strouhal number for flow past a circular cylinder. If the flow is steady, only drag and lift coefficients are computed. Call signature: draglift(simulations, d_cylinder=0.1, u_0=1.0, flow_dir='y', t_start=-1, sortby='dx'): Keyword arguments: *simulations* array of simulation names to be included in the computations *d_cylinder*: diameter of the cylinder *u_0* velocity at the inlet *flow_dir*: direction of the flow *t_start* time to start the drag computations from should be where the a steady vortex shedding has developed *sortby* property to sort the arrays by typical choices for parametric studies are grid size, length of domain, etc. Returns object with information: sim-name, drag (mean), lift (rms), and st (non-dim shedding frequency) .. py:class:: Draglift Bases: :py:obj:`object` Grid -- holds pencil code time grid data. Fill members with default values. .. py:attribute:: name :value: '' .. py:attribute:: drag :value: 0 .. py:attribute:: lift :value: 0 .. py:attribute:: st :value: 0 .. py:method:: set_name(simulation) .. py:method:: compute(simulation, d_cylinder=0.1, u_0=1.0, flow_dir='y', t_start=-1) Compute the drag coefficienc, rms drag, lift coefficient, rms lift and Strouhal number of a given time series Requires time series T as input, on the form given in the pencil code where T.t, T.c_dragx and T.c_dragy are avaliable quantities Cylinder diameter, fluid mean velocity , flow direction and start time of drag computations velocity should also be given as input. .. py:function:: find_extrema(series, maxmin) Input: Data series, extrama of intereset (max or min) Use scipy to find this scipy.signal.argrelextrema(array type of np,comparison operator eg np.greater or np.less) .. py:function:: grav_profile(grav, x, y, z, par=None) .. py:function:: tensors_sph(*args, **kwargs) .. py:class:: Tensors Bases: :py:obj:`object` Tensors -- Holds the calculated z-averaged tensors and code coefficients Fill members with default values. .. py:attribute:: t .. py:method:: calc(aver=[], datatopdir='.', lskip_zeros=False, proc=-1, rank=0, rmfzeros=1, rmbzeros=1, iy=None, l_correction=False, t_correction=0.0, dim=None, timereducer=None, trargs=[], tindex=(0, None), imask=None) object returns time dependent meridional tensors from Averages object aver.z. u, acoef and bcoef and aver.t For long DNS runs the 'zaverages.dat' file can be very large so MPI may be required and the data is loaded by processor as default. lskip_zeros=True identifies the resetting of the testfield and rmbzeros and rmfzeros number to exclude before and following By default none are removed. iy is the index array that is computed in this MPI process, which may be a subset of the array on this processor l_correction=True permits the pencil coefficients computed prior to the Pencil Code correction implemented after time=t_correction to be rescaled accordingly to match the new formulation. trargs contain optional arguments for the time treatments: mean, smoothing, etc. tindex is set to limit the range of the iterations loaded from Averages in zaverages.dat The index imask, excluding the resets, can be specified to ensure all processes use the same mask .. py:function:: dot(a, b) dot(a, b) Take dot product of two pencil-code vectors a and b. :param a: Pencil-code vectors with shape [3, mz, my, mx]. :type a: ndarrays :param b: Pencil-code vectors with shape [3, mz, my, mx]. :type b: ndarrays .. py:function:: dot2(a) dot2(a) Take dot product of a pencil-code vector with itself. :param a: Pencil-code vector with shape [3, mz, my, mx]. :type a: ndarray .. py:function:: cross(a, b) cross(a, b) Take cross of two pencil-code vectors a and b. :param a: Pencil-code vectors with shape [3, mz, my, mx]. :type a: ndarrays :param b: Pencil-code vectors with shape [3, mz, my, mx]. :type b: ndarrays .. py:function:: div(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system='cartesian', grid=None) div(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system="cartesian", grid=None) Take divergence of pencil code vector array f in various coordinate systems. :param f: Pencil code vector array f. :type f: ndarray :param grid: Pencil grid object. See pc.read.grid(). :type grid: pencil.read.grids.Grid :param coordinate_system: Coordinate system under which to take the divergence. Takes 'cartesian', 'cylindrical' and 'spherical'. :type coordinate_system: string :param Deprecated parameters (only for backwards compatibility): :param --------------------------------------------------------: :param dx: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dx: floats :param dy: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dy: floats :param dz: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dz: floats :param x: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type x: ndarrays :param y: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type y: ndarrays .. py:function:: curl(f, dx=None, dy=None, dz=None, x=None, y=None, run2D=False, coordinate_system='cartesian', grid=None) curl(f, dx=None, dy=None, dz=None, x=None, y=None, run2D=False, coordinate_system="cartesian", grid=None) Take the curl of a pencil code vector array f in various coordinate systems. :param f: Pencil code scalar array f. :type f: ndarray :param grid: Pencil grid object. See pc.read.grid(). :type grid: pencil.read.grids.Grid :param run2D: Deals with pure 2-D snapshots. !Only for Cartesian grids at the moment! Requires grid!=None. :type run2D: bool :param coordinate_system: Coordinate system under which to take the divergence. Takes 'cartesian', 'cylindrical' and 'spherical'. !Does not work for 2d runs yet! :type coordinate_system: string :param Deprecated parameters (only for backwards compatibility): :param --------------------------------------------------------: :param dx: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dx: floats :param dy: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dy: floats :param dz: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dz: floats :param x: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type x: ndarrays :param y: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type y: ndarrays .. py:function:: curl2(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system='cartesian', grid=None) curl2(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system="cartesian", grid=None) Take the double curl of a pencil code vector array f. :param f: Pencil code vector array f. :type f: ndarray :param grid: Pencil grid object. See pc.read.grid(). :type grid: pencil.read.grids.Grid :param coordinate_system: Coordinate system under which to take the divergence. Takes 'cartesian', 'cylindrical' and 'spherical'. :type coordinate_system: string :param Deprecated parameters (only for backwards compatibility): :param --------------------------------------------------------: :param dx: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dx: floats :param dy: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dy: floats :param dz: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dz: floats :param x: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type x: ndarrays :param y: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type y: ndarrays .. py:function:: grad(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system='cartesian', grid=None) grad(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system="cartesian", grid=None) Take the gradient of a pencil code scalar array f in various coordinate systems. :param f: Pencil code scalar array f. :type f: ndarray :param grid: Pencil grid object. See pc.read.grid(). :type grid: pencil.read.grids.Grid :param coordinate_system: Coordinate system under which to take the divergence. Takes 'cartesian', 'cylindrical' and 'spherical'. :type coordinate_system: string :param Deprecated parameters (only for backwards compatibility): :param --------------------------------------------------------: :param dx: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dx: floats :param dy: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dy: floats :param dz: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dz: floats :param x: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type x: ndarrays :param y: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type y: ndarrays .. py:function:: del2(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system='cartesian', grid=None) del2(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system="cartesian", grid=None) Calculate del2, the Laplacian of a scalar field f. :param f: Pencil code vector array f. :type f: ndarray :param grid: Pencil grid object. See pc.read.grid(). :type grid: pencil.read.grids.Grid :param coordinate_system: Coordinate system under which to take the divergence. Takes 'cartesian', 'cylindrical' and 'spherical'. :type coordinate_system: string :param Deprecated parameters (only for backwards compatibility): :param --------------------------------------------------------: :param dx: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dx: floats :param dy: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dy: floats :param dz: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dz: floats :param x: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type x: ndarrays :param y: Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None :type y: ndarrays .. py:function:: del6(f, dx=None, dy=None, dz=None, grid=None) del6(f, dx=None, dy=None, dz=None, grid=None) Calculate del6 (defined here as d^6/dx^6 + d^6/dy^6 + d^6/dz^6, rather than del2^3) of a scalar f for hyperdiffusion. :param f: Pencil code scalar array f. :type f: ndarray :param grid: Pencil grid object. See pc.read.grid(). :type grid: pencil.read.grids.Grid :param Deprecated parameters (only for backwards compatibility): :param --------------------------------------------------------: :param dx: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dx: floats :param dy: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dy: floats :param dz: Grid spacing in the three dimensions. These will not have any effect if grid!=None :type dz: floats .. py:function:: fluid_reynolds(uu, param, grid, lnrho=list(), shock=list(), nghost=3, lmix=True, quiet=True) Computes the fluid Reynolds number from the advective and effective viscous expressions in the momentum equation. call signature: fluid_reynolds(uu, ivisc, grid, rho=None, shock=None, nghost=3) :keyword \*uu*: The velocity field [3,mz,my,mx] from the simulation data :keyword \*param*: The Param simulation object with viscosity data information :keyword \*grid*: The Grid simulation object :keyword \*lnrho*: The log density field if it is non-uniform :keyword \*shock*: The shock variable if shock viscosity is applied :keyword \*nghost*: The number of ghost zones appropriate to the order of accuracy :keyword \*lmix*: Option not to include hyper values when Laplacian values present .. py:function:: magnetic_reynolds(uu, param, grid, aa=list(), bb=list(), jj=list(), nghost=3, lmix=True, quiet=True) Computes the magnetic Reynolds number from the advective and effective resistive expressions in the induction equation. call signature: magnetic_reynolds(uu, param, grid, aa=None, bb=None, jj=None, nghost=3): :keyword \*uu*: The velocity field [3,mz,my,mx] from the simulation data :keyword \*param*: The Param simulation object with resistivity data information :keyword \*grid*: The Grid simulation object :keyword \*aa*: The vector potential if bb is not present or hyper diffusion :keyword \*bb*: The magnetic field :keyword \*jj*: The current density field :keyword \*nghost*: The number of ghost zones appropriate to the order of accuracy :keyword \*lmix*: Option not to include hyper values when Laplacian values present .. py:function:: sod(*args, **kwargs) Solve Sod shock tube for given parameters and resolution Signature: calc_shocktube(xarr, time, par=list(), lreference=False, DEBUG=False, lplot=False, itplot=0, magic=['ee','tt','Ms']) :param \*xarr*: :type \*xarr*: Coordinate vector (the initial discontinuity is always at x=0) :param \*time*: :type \*time*: Time after membrane snapped. :param \*par*: :type \*par*: Param object :param \*lreference*: :type \*lreference*: Use default parameters for Sod's reference problem. :param \*DEBUG*: :type \*DEBUG*: Flag to switch on output :param \*lplot*: :type \*lplot*: Plot first snapshot profiles :param \*itplot*: :type \*itplot*: Iteration index of snaphot to plot :param \*magic*: :type \*magic*: Optional profiles to include in Sod object. :rtype: Class containing coordinates and shock profiles .. rubric:: Notes Adapted from IDL script shocktube.pro Analytical solution of the shocktube problem Initial state (pressure jump) ------------------+------------------ pl | pr ------------------+------------------ 0 Evolved state: membrane has snapped, four domains ----------------------------------------- 1 (l) | 2 | 3 | 4 | 5 (r) --------------+------+-------+----+------ x x x x 1 2 3 4 Velocity ************* * * * -**************--------------------***--- Nomenclature pl, pr - pressure far left/right from the shock structure p2 - pressure in the expansion wave p3=p4 - no pressure jump across the contact discontinuity p: pressure, T: temperature, rho: density, u: flow velocity, cs: sound speed gamma (adiabatic index) is assumed to be constant The borders between the different domains are x1, x2, x3, x4 and move at velocities ul-csl, u4-c4, u2, u3, u4, respectively. Warning: This works so far only in the case ul=ur=0 .. rubric:: Examples >>> shock = pc.calc.sod() >>> shock.keys() t x rho ux pp ee tt Ms .. py:function:: kernel_smooth(sim_path, src, dst, magic=['meanuu'], par=[], comm=None, gd=[], grp_overwrite=False, overwrite=False, rank=0, size=1, nghost=3, kernel=1.0, status='a', chunksize=1000.0, dtype=np.float64, quiet=True, nmin=32, typ='piecewise', mode=list()) .. py:function:: gauss_3Dsmooth(arr, sigma=1.0, typ='all', quiet=True, mode=list())