/usr/lib/python2.7/dist-packages/PySPH-1.0a4.dev0-py2.7-linux-x86_64.egg/pysph/examples/rayleigh_taylor.py is in python-pysph 0~20160514.git91867dc-4build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 | """Rayleigh-Taylor instability problem. (16 hours)
"""
# numpy
import numpy as np
# PySPH imports
from pysph.base.utils import get_particle_array
from pysph.solver.solver import Solver
from pysph.solver.application import Application
from pysph.sph.scheme import TVFScheme
# domain and reference values
gy = -1.0
Lx = 1.0; Ly = 2.0
Re = 420; Vmax = np.sqrt(0.5*Ly*abs(gy))
nu = Vmax*Ly/Re
# density for the two phases
rho1 = 1.8; rho2 = 1.0
# speed of sound and reference pressure
Fr = 0.01
c0 = Vmax/Fr
p1 = c0**2 * rho1
p2 = c0**2 * rho2
# Numerical setup
nx = 50; dx = Lx/nx
ghost_extent = 5 * dx
hdx = 1.2
# adaptive time steps
h0 = hdx * dx
dt_cfl = 0.25 * h0/( c0 + Vmax )
dt_viscous = 0.125 * h0**2/nu
dt_force = 0.25 * np.sqrt(h0/abs(gy))
tf = 25
dt = 0.5 * min(dt_cfl, dt_viscous, dt_force)
class RayleighTaylor(Application):
def create_particles(self):
# create all the particles
_x = np.arange( -ghost_extent - dx/2, Lx + ghost_extent + dx/2, dx )
_y = np.arange( -ghost_extent - dx/2, Ly + ghost_extent + dx/2, dx )
x, y = np.meshgrid(_x, _y); x = x.ravel(); y = y.ravel()
# sort out the fluid and the solid
indices = []
for i in range(x.size):
if ( (x[i] > 0.0) and (x[i] < Lx) ):
if ( (y[i] > 0.0) and (y[i] < Ly) ):
indices.append(i)
# create the arrays
solid = get_particle_array(name='solid', x=x, y=y)
# remove the fluid particles from the solid
fluid = solid.extract_particles(indices); fluid.set_name('fluid')
solid.remove_particles(indices)
# sort out the two fluid phases
indices = []
for i in range(fluid.get_number_of_particles()):
if fluid.y[i] > 1 - 0.15*np.sin(2*np.pi*fluid.x[i]):
indices.append(i)
fluid1 = fluid.extract_particles(indices); fluid1.set_name('fluid1')
fluid2 = fluid
fluid2.set_name('fluid2')
fluid2.remove_particles(indices)
fluid1.rho[:] = rho1
fluid2.rho[:] = rho2
print("Rayleigh Taylor Instability problem :: Re = %d, nfluid = %d, nsolid=%d, dt = %g"%(
Re, fluid1.get_number_of_particles() + fluid2.get_number_of_particles(),
solid.get_number_of_particles(), dt))
# add requisite properties to the arrays:
self.scheme.setup_properties([fluid1, fluid2, solid])
# setup the particle properties
volume = dx * dx
# mass is set to get the reference density of each phase
fluid1.m[:] = volume * rho1
fluid2.m[:] = volume * rho2
# volume is set as dx^2
fluid1.V[:] = 1./volume
fluid2.V[:] = 1./volume
solid.V[:] = 1./volume
# smoothing lengths
fluid1.h[:] = hdx * dx
fluid2.h[:] = hdx * dx
solid.h[:] = hdx * dx
# return the arrays
return [fluid1, fluid2, solid]
def create_scheme(self):
s = TVFScheme(
['fluid1', 'fluid2'], ['solid'], dim=2, rho0=rho1, c0=c0, nu=nu,
p0=p1, pb=p1, h0=dx*hdx, gy=gy
)
s.configure_solver(tf=tf, dt=dt, pfreq=500)
return s
def create_equations(self):
# This is an ugly hack to support different densities for fluids.
# What we should really do is set rho0 as a fluid constant and rewrite
# the equations to use that, then once the fluid properties are
# defined, this will just work.
equations = super(RayleighTaylor, self).create_equations()
from pysph.sph.equation import Group
def process_term(x):
if hasattr(x, 'rho0'):
if x.dest == 'fluid1' or x.sources == ['fluid1']:
x.rho0 = rho1
elif x.dest == 'fluid2' or x.sources == ['fluid2']:
x.rho0 = rho2
if hasattr(x, 'p0'):
if x.dest == 'fluid1':
x.p0 = p1
elif x.dest == 'fluid2':
x.p0 = p2
for eq in equations:
if isinstance(eq, Group):
for eq in eq.equations:
process_term(eq)
else:
process_term(eq)
return equations
if __name__ == '__main__':
app = RayleighTaylor()
app.run()
|