python-pysph 0~20160514.git91867dc-4build1 / usr / lib / python2.7 / dist-packages / PySPH-1.0a4.dev0-py2.7-linux-x86_64.egg / pysph / examples / rayleigh_taylor.py

This file is indexed.

/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()

APT Browse - Built by Thomas Orozco.