pencil.calc

Math functions and further calculations.

Submodules

Attributes

pc_print

Classes

Stream

Contains the methods and results for the streamline tracing for a field on a grid.

Draglift

Grid -- holds pencil code time grid data.

Tensors

Tensors -- Holds the calculated z-averaged tensors and code

Functions

part_to_grid(xp, yp[, zp, quantity, Nbins, sim, ...])

Bins quantity based on position data xp and yp to 1024^2 bins like a histrogram.

fill_gaps_in_grid(array[, key, steps, order, DEBUG])

Interpolates nans by nearest neighbor method to fill gaps in arrays.

twonorm_accuracy(simulations[, field, strip, ...])

Assessment of accuracy of simulation:

order_accuracy([simulations, nstrips, twonorm_arr, ...])

Compute an estimate of the order of accuracy, using two-norms where

twonorm_array(simulations, nstrips[, field, var_file, ...])

Compute twonorm_accuracy along the selected direction, for a number of 'strips'.

twonorm(u1, u2, dx)

Compute the two-norm error of two spatial vectors u1 and u2.

maxerror(u1, u2)

Find the larges error between values from u1 and u2 (for the

twonorm_accuracy2D(simulations[, field, var_file, ...])

Assessment of accuracy of simulation:

twonorm_accuracy1D(simulations[, field, strip, ...])

Assessment of accuracy of simulation:

draglift(simulations[, sortby])

Compute mean drag coefficient, rms lift coefficient, and Strouhal number

find_extrema(series, maxmin)

Input: Data series, extrama of intereset (max or min)

grav_profile(grav, x, y, z[, par])

tensors_sph(*args, **kwargs)

dot(a, b)

dot(a, b)

dot2(a)

dot2(a)

cross(a, b)

cross(a, b)

div(f[, dx, dy, dz, x, y, coordinate_system, grid])

div(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system="cartesian", grid=None)

curl(f[, dx, dy, dz, x, y, run2D, coordinate_system, grid])

curl(f, dx=None, dy=None, dz=None, x=None, y=None, run2D=False, coordinate_system="cartesian", grid=None)

curl2(f[, dx, dy, dz, x, y, coordinate_system, grid])

curl2(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system="cartesian", grid=None)

grad(f[, dx, dy, dz, x, y, coordinate_system, grid])

grad(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system="cartesian", grid=None)

del2(f[, dx, dy, dz, x, y, coordinate_system, grid])

del2(f, dx=None, dy=None, dz=None, x=None, y=None, coordinate_system="cartesian", grid=None)

del6(f[, dx, dy, dz, grid])

del6(f, dx=None, dy=None, dz=None, grid=None)

fluid_reynolds(uu, param, grid[, lnrho, shock, ...])

Computes the fluid Reynolds number from the advective and effective

magnetic_reynolds(uu, param, grid[, aa, bb, jj, ...])

Computes the magnetic Reynolds number from the advective and effective

sod(*args, **kwargs)

Solve Sod shock tube for given parameters and resolution

kernel_smooth(sim_path, src, dst[, magic, par, comm, ...])

gauss_3Dsmooth(arr[, sigma, typ, quiet, mode])

Package Contents

pencil.calc.pc_print
class pencil.calc.Stream(field, params, xx=(0, 0, 0), time=(0, 1), metric=None, splines=None)

Bases: 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)])

params
xx = (0, 0, 0)
time = (0, 1)
section_l
total_l
iterations
section_dh
total_h = 1
pencil.calc.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.

Parameters:

  • xp (-) – array of x and y positions

  • yp – array of x and y positions

  • zp (-) – specify if 3D run, set False to have zp == 0

  • quantity (-) – array of same shape as xp and yp, but with quantity to bin, set it False to count number of occurrences/histrogram2d

  • Nbins (-) – number of histrogram bins for each direction. if 2d only the first two entries in Nbins are used

  • sim (-) – to extract extent from by loading ghostzone-free grid

  • 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

  • cells (- fill_gaps interpolate empty grid)

Returns: arr, xgrid, ygrid

  • arr: 2d array with binned values

  • x-/ygrid: linspace of used x/y grid

  • zgrid: if zp != False

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

pencil.calc.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!

Parameters:

  • array (-) – array with nans inside

  • key (-) – indicates array entries that will be interpolated

  • method (-) – as specified in scipy.interpolate.griddata, i.e. use ‘linear’ for 3D data

Returns:

array without nans inside.

Example

[[0, 1, 0], [[0, 1, 0,],

[1, nan, 1], -interpol-> [1, 2, 1], [0, 1, 0]] [0, 1, 0]]

pencil.calc.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 sumlimits_{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

pencil.calc.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

pencil.calc.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

pencil.calc.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.

pencil.calc.maxerror(u1, u2)

Find the larges error between values from u1 and u2 (for the same indices).

pencil.calc.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 sumlimits_{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

pencil.calc.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 sumlimits_{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

pencil.calc.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)

class pencil.calc.Draglift

Bases: object

Grid – holds pencil code time grid data.

Fill members with default values.

name = ''
drag = 0
lift = 0
st = 0
set_name(simulation)
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.

pencil.calc.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)

pencil.calc.grav_profile(grav, x, y, z, par=None)
pencil.calc.tensors_sph(*args, **kwargs)
class pencil.calc.Tensors

Bases: object

Tensors – Holds the calculated z-averaged tensors and code coefficients

Fill members with default values.

t
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

pencil.calc.dot(a, b)

dot(a, b)

Take dot product of two pencil-code vectors a and b.

Parameters:

  • a (ndarrays) – Pencil-code vectors with shape [3, mz, my, mx].

  • b (ndarrays) – Pencil-code vectors with shape [3, mz, my, mx].

pencil.calc.dot2(a)

dot2(a)

Take dot product of a pencil-code vector with itself.

Parameters:

a (ndarray) – Pencil-code vector with shape [3, mz, my, mx].

pencil.calc.cross(a, b)

cross(a, b)

Take cross of two pencil-code vectors a and b.

Parameters:

  • a (ndarrays) – Pencil-code vectors with shape [3, mz, my, mx].

  • b (ndarrays) – Pencil-code vectors with shape [3, mz, my, mx].

pencil.calc.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.

Parameters:

  • f (ndarray) – Pencil code vector array f.

  • grid (pencil.read.grids.Grid) – Pencil grid object. See pc.read.grid().

  • coordinate_system (string) – Coordinate system under which to take the divergence. Takes ‘cartesian’, ‘cylindrical’ and ‘spherical’.

  • compatibility) (Deprecated parameters (only for backwards)

  • --------------------------------------------------------

  • dx (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dy (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dz (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • x (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

  • y (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

pencil.calc.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.

Parameters:

  • f (ndarray) – Pencil code scalar array f.

  • grid (pencil.read.grids.Grid) – Pencil grid object. See pc.read.grid().

  • run2D (bool) – Deals with pure 2-D snapshots. !Only for Cartesian grids at the moment! Requires grid!=None.

  • coordinate_system (string) – Coordinate system under which to take the divergence. Takes ‘cartesian’, ‘cylindrical’ and ‘spherical’. !Does not work for 2d runs yet!

  • compatibility) (Deprecated parameters (only for backwards)

  • --------------------------------------------------------

  • dx (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dy (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dz (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • x (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

  • y (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

pencil.calc.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.

Parameters:

  • f (ndarray) – Pencil code vector array f.

  • grid (pencil.read.grids.Grid) – Pencil grid object. See pc.read.grid().

  • coordinate_system (string) – Coordinate system under which to take the divergence. Takes ‘cartesian’, ‘cylindrical’ and ‘spherical’.

  • compatibility) (Deprecated parameters (only for backwards)

  • --------------------------------------------------------

  • dx (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dy (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dz (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • x (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

  • y (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

pencil.calc.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.

Parameters:

  • f (ndarray) – Pencil code scalar array f.

  • grid (pencil.read.grids.Grid) – Pencil grid object. See pc.read.grid().

  • coordinate_system (string) – Coordinate system under which to take the divergence. Takes ‘cartesian’, ‘cylindrical’ and ‘spherical’.

  • compatibility) (Deprecated parameters (only for backwards)

  • --------------------------------------------------------

  • dx (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dy (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dz (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • x (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

  • y (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

pencil.calc.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.

Parameters:

  • f (ndarray) – Pencil code vector array f.

  • grid (pencil.read.grids.Grid) – Pencil grid object. See pc.read.grid().

  • coordinate_system (string) – Coordinate system under which to take the divergence. Takes ‘cartesian’, ‘cylindrical’ and ‘spherical’.

  • compatibility) (Deprecated parameters (only for backwards)

  • --------------------------------------------------------

  • dx (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dy (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dz (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • x (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

  • y (ndarrays) – Radial (x) and polar (y) coordinates, 1d arrays. These will not have any effect if grid!=None

pencil.calc.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.

Parameters:

  • f (ndarray) – Pencil code scalar array f.

  • grid (pencil.read.grids.Grid) – Pencil grid object. See pc.read.grid().

  • compatibility) (Deprecated parameters (only for backwards)

  • --------------------------------------------------------

  • dx (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dy (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

  • dz (floats) – Grid spacing in the three dimensions. These will not have any effect if grid!=None

pencil.calc.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 Arguments:

  • *uu* – The velocity field [3,mz,my,mx] from the simulation data

  • *param* – The Param simulation object with viscosity data information

  • *grid* – The Grid simulation object

  • *lnrho* – The log density field if it is non-uniform

  • *shock* – The shock variable if shock viscosity is applied

  • *nghost* – The number of ghost zones appropriate to the order of accuracy

  • *lmix* – Option not to include hyper values when Laplacian values present

pencil.calc.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 Arguments:

  • *uu* – The velocity field [3,mz,my,mx] from the simulation data

  • *param* – The Param simulation object with resistivity data information

  • *grid* – The Grid simulation object

  • *aa* – The vector potential if bb is not present or hyper diffusion

  • *bb* – The magnetic field

  • *jj* – The current density field

  • *nghost* – The number of ghost zones appropriate to the order of accuracy

  • *lmix* – Option not to include hyper values when Laplacian values present

pencil.calc.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’])

Parameters:

  • *xarr* (Coordinate vector (the initial discontinuity is always at x=0))

  • *time* (Time after membrane snapped.)

  • *par* (Param object)

  • *lreference* (Use default parameters for Sod's reference problem.)

  • *DEBUG* (Flag to switch on output)

  • *lplot* (Plot first snapshot profiles)

  • *itplot* (Iteration index of snaphot to plot)

  • *magic* (Optional profiles to include in Sod object.)

Return type:

Class containing coordinates and shock profiles

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

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

Examples

>>> shock = pc.calc.sod()
>>> shock.keys()
t
x
rho
ux
pp
ee
tt
Ms
pencil.calc.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())
pencil.calc.gauss_3Dsmooth(arr, sigma=1.0, typ='all', quiet=True, mode=list())