def print_strengths(sid, fluid=False, c=0.0):
    sim = sphere.sim(sid, fluid=fluid)

    sim.readfirst(verbose=False)
    sim.visualize('shear')

    friction = sim.tau[1:]/sim.sigma_eff[1:]
    tau_peak = numpy.max(friction)
    tau_ultimate = numpy.average(friction[-500:-1])

    if fluid:
        print('%.2f \t %.2f \t %.2f' % (c, tau_peak, tau_ultimate))
    else:
        print('dry \t %.2f \t %.2f' % (tau_peak, tau_ultimate))

    return friction

#!/usr/bin/env python
from pytestutils import *
import sphere

#### Input/output tests ####
print("### Fluid input/output tests ###")

# Generate data in python
orig = sphere.sim(np=100, sid="test-initgrid-fluid", fluid=True)
orig.generateRadii(histogram=False, radius_mean=1.0)
orig.defaultParams()
orig.initRandomGridPos()
orig.initFluid()
orig.initTemporal(current=0.0, total=0.0)
orig.time_total=2.0*orig.time_dt
orig.time_file_dt = orig.time_dt
orig.writebin(verbose=False)

# Test Python IO routines
py = sphere.sim(fluid=True)
py.readbin("../input/" + orig.sid + ".bin", verbose=False)
compare(orig, py, "Python IO:")

# Test C++ IO routines
orig.run()
#orig.run(dry=True)
#orig.run(verbose=True, hideinputfile=False, cudamemcheck=True)
cpp = sphere.sim(fluid=True)
cpp.readbin("../output/" + orig.sid + ".output00000.bin", verbose=False)
compare(orig, cpp, "C++ IO:   ")

if smoothed_results:
    friction_smooth = [[], [], [], []]
dilation = [[], [], [], []]
p_min = [[], [], [], []]
p_mean = [[], [], [], []]
p_max = [[], [], [], []]
f_n_mean = [[], [], [], []]
f_n_max  = [[], [], [], []]
v_f_z_mean  = [[], [], [], []]

fluid=True

# dry shear
#sid = 'shear-sigma0=' + sys.argv[1] + '-hw'
sid = 'halfshear-sigma0=' + sys.argv[1] + '-shear'
sim = sphere.sim(sid)
sim.readlast(verbose=False)
sim.visualize('shear')
shear_strain[0] = sim.shear_strain
#shear_strain[0] = numpy.arange(sim.status()+1)
friction[0] = sim.tau/sim.sigma_eff
if smoothed_results:
    friction_smooth[0] = smooth(friction[0], smooth_window)
dilation[0] = sim.dilation

if contact_forces:
    f_n_mean[0] = numpy.zeros_like(shear_strain[0])
    f_n_max[0]  = numpy.zeros_like(shear_strain[0])
    for i in numpy.arange(sim.status()):
        sim.readstep(i, verbose=False)
        sim.findNormalForces()

#!/usr/bin/env python
import sphere
import numpy


for dp in [1.0e3, 2.0e3]:
    # Read initial configuration
    sim = sphere.sim('diffusivity-relax')
    sim.readlast()

    sim.sid = 'permeability-dp=' + str(dp)
    sim.cleanup()

    sim.g[2] = 0.0
    sim.nw[0] = 0
    sim.initGrid()
    sim.zeroKinematics()
    sim.initFluid(mu = 17.87e-4, p = 1.0e5, hydrostatic=True)

    # Initialize to linear hydraulic gradient
    p_bottom = 10.0
    p_top = p_bottom + dp
    dz = sim.L[2]/sim.num[2]
    for iz in range(sim.num[2]-1):
        #z = dz*iz + 0.5*dz # cell-center z-coordinate
        z = dz*iz
        sim.p_f[:,:,iz] = p_bottom + dp/sim.L[2] * z

    sim.setFluidTopFixedPressure()
    sim.setFluidBottomFixedPressure()
    sim.p_f[:,:,-1] = p_top

#!/usr/bin/env python
import sphere
import numpy

cube = sphere.sim('cube-init')
cube.readlast()
cube.adjustUpperWall(z_adjust=1.0)

# Fill out grid with cubic packages
grid = numpy.array((
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2],
        [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
        [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
        [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
        [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
        [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
        [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
        [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
        [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
        [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
        ))

# CUDA device to use
device = 0

# Number of particles
np = 1e4

# Common simulation id
sim_id = "jp-long-shear-soft"

# Deviatoric stress [Pa]
devslist = [100e3]

### INITIALIZATION ###

# New class
init = sphere.sim(np=np, nd=3, nw=0, sid=sim_id + "-init")

# Save radii
init.generateRadii(mean=0.01)

# Use default params
init.defaultParams(k_n=1.16e7, k_t=1.16e7, gamma_n=100.0, mu_s=0.6, mu_d=0.6)

# Add gravity
init.g[2] = -9.81

# Periodic x and y boundaries
init.periodicBoundariesXY()

# Initialize positions in random grid (also sets world size)
hcells = np**(1.0/3.0)

    i = 0

    for sid in sids:
        if os.path.isfile('../output/' + sid + '.status.dat'):
            pc = PermeabilityCalc(sid, plot_evolution=False,
                    print_results=False, verbose=False)
            K[c][i] = pc.conductivity()
            pc.findPressureGradient()
            pc.findCrossSectionalFlux()
            dpdz[c][i] = pc.dPdL[2]
            Q[c][i] = pc.Q[2]
            pc.findMeanPorosity()
            #pc.plotEvolution()
            phi_bar[c][i] = pc.phi_bar

            sim = sphere.sim(sid, fluid=True)
            sim.readlast(verbose=False)
            Re[c][i] = numpy.mean(sim.ReynoldsNumber())

            #sim.writeVTKall()

            # find magnitude of fluid pressure force and total interaction force
            '''
            fp_magn = numpy.empty(sim.np)
            fsum_magn = numpy.empty(sim.np)
            for i in numpy.arange(sim.np):
                fp_magn[i] = sim.f_p[i,:].dot(sim.f_p[i,:])
                fsum_magn[i] = sim.f_sum[i,:].dot(sim.f_sum[i,:])

            fp_fsum[c][i] = numpy.mean(fp_magn/fsum_magn)
            # interaction forces not written in these old output files!

#!/usr/bin/env python
import sphere

sim = sphere.sim('init2', np=10000)
sim.generateRadii(psd='uni', radius_mean=0.02, radius_variance=0.01)
sim.initRandomGridPos([12, 12, 1000])
sim.initTemporal(10.0, file_dt=0.05, epsilon=0.07)
sim.gamma_n[0] = 1000.0
sim.gamma_wn[0] = 1000.0
sim.periodicBoundariesXY()
sim.g[2] = -9.81
sim.run()
sim.writeVTKall()

#!/usr/bin/env python
from pytestutils import *
import sphere

#### Input/output tests ####
print("### Memory tests ###")

# Generate data in python
orig = sphere.sim(np = 100, nw = 1, sid = "test-initgrid")
orig.generateRadii(histogram = False)
orig.defaultParams()
orig.initRandomGridPos()
orig.initTemporal(current = 0.0, total = 0.0)
orig.time_total = 2.0*orig.time_dt;
orig.time_file_dt = orig.time_dt;

# Test C++ routines
print("Valgrind: C++ routines")
orig.run(verbose=False, hideinputfile=True, valgrind=True)


# Test CUDA routines
print("cuda-memcheck: CUDA routines")
orig.run(verbose=False, hideinputfile=True, cudamemcheck=True)

# Remove temporary files
cleanup(orig)

# cohesion, 1 denotes a wet simulation with capillary cohesion.
# GRAVITY toggles gravitational acceleration. Without it, the particles are
# placed in the middle of a volume. With it enabled, the particles are put on
# top of a flat wall.

import sphere

# import numpy
import sys

device = int(sys.argv[1])
cohesion = int(sys.argv[2])
gravity = int(sys.argv[3])

# Create packing
sim = sphere.sim("cap-cohesion=" + str(cohesion) + "-init-grav=" + str(gravity), np=2000)
# sim.mu_s[0] = 0.0
# sim.mu_d[0] = 0.0
# sim.k_n[0] = 1.0e7
# sim.k_t[0] = 1.0e7
sim.generateRadii(psd="uni", mean=1.0e-3, variance=1.0e-4)
sim.contactModel(1)
sim.initRandomGridPos(gridnum=[24, 24, 10000], padding=1.4)
sim.defaultParams(gamma_t=1.0e3, capillaryCohesion=1)
sim.initTemporal(5.0, file_dt=0.01, epsilon=0.07)
# I = numpy.nonzero(sim.x[:,2] < sim.L[2]*0.5)
# sim.vel[I[0], 2] =  0.01  # add a instability seeding perturbation
# I = numpy.nonzero(sim.x[:,2] > sim.L[2]*0.5)
# sim.vel[I[0], 2] = -0.01  # add a instability seeding perturbation
if gravity == 1:
    sim.g[2] = -10.0

#!/usr/bin/env python
import sphere
import numpy
import subprocess

sid = 'shear-sigma0=20000.0-hw'
imgformat = 'png'
sim = sphere.sim(sid, fluid=False)
subprocess.call('mkdir -p ' + sid + '-fc', shell=True)

d = 0
for i in numpy.arange(210,300):
#for i in numpy.arange(240,260):
    print("File: %d, output: %s-fc/%05d.png" % (i, sid, d))
    sim.readstep(i, verbose=False)
    #sim.forcechains(lc = 3.0e1, uc=1.0e2)
    sim.forcechains(lc = 30.0, uc=1.0e2)
    subprocess.call('mv shear-sigma0=20000-0-hw-fc.' + imgformat \
            + ' ' + sid + '-fc/%05d.png' % (d), shell=True)
    d += 1

subprocess.call('cd ' + sid + '-fc && sh ./make_video.sh', shell=True)

#!/usr/bin/env python
import sphere
import numpy
import sys

# launch with:
# $ python consolidation-curves <DEVICE> <C_PHI> <C_GRAD_P>

# Unique simulation parameters
device = int(sys.argv[1])
c_phi = float(sys.argv[2])
c_grad_p = float(sys.argv[3])

sim = sphere.sim('cons-sigma0=' + str(5.0e3) + '-c_phi=' + \
            str(c_phi) + '-c_grad_p=' + str(c_grad_p), fluid=True)
sim.readlast()

sigma0 = 10.0e3
sim.sid = 'cons-sigma0=' + str(sigma0) + '-c_phi=' + str(c_phi) + \
        '-c_grad_p=' + str(c_grad_p) + '-tall'
print('\n###### ' + sim.sid + ' ######')

# Checkerboard colors
sim.checkerboardColors()
sim.cleanup()
#sim.adjustUpperWall()
sim.zeroKinematics()

#sim.consolidate(normal_stress = sigma0)
sim.w_sigma0[0] = sigma0

# Fluid density
rho_f = 1000.0

# Gravitational acceleration
g = 10.0

# Number of particles
np = 1e4


### INITIALIZATION ###

# New class
# init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
init = sphere.sim(np=np, nd=3, nw=0, sid="creep1-init")

# Uniform radii from 0.8 cm to 1.2 cm
init.generateRadii(psd="uni", mean=0.005, variance=0.001)

# Use default params
init.defaultParams(gamma_n=100.0, mu_s=0.6, mu_d=0.6)
init.setYoungsModulus(1e8)

# Add gravity
init.g[2] = -g

# Periodic x and y boundaries
init.periodicBoundariesXY()

# Initialize positions in random grid (also sets world size)

#!/usr/bin/env python
from pytestutils import *

import sphere
import sys
import numpy
import matplotlib.pyplot as plt

print('### CFD tests - Dirichlet/Neumann BCs ###')

print('''# Neumann bottom, Dirichlet top BC.
# No gravity, no pressure gradients => no flow''')
orig = sphere.sim("neumann", fluid = True)
cleanup(orig)
orig.defaultParams(mu_s = 0.4, mu_d = 0.4)
orig.defineWorldBoundaries([0.4, 0.4, 1], dx = 0.1)
orig.initFluid(mu = 8.9e-4)
#orig.initFluid(mu = 0.0)
orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
orig.c_v[0] = 0.1
#orig.c_phi[0] = 0.1
py = sphere.sim(sid = orig.sid, fluid = True)
orig.bc_bot[0] = 1      # No-flow BC at bottom (Neumann)
#orig.run(dry=True)
orig.run(verbose=False)
#orig.run(device=2)
#orig.writeVTKall()
py.readlast(verbose = False)
zeros = numpy.zeros((orig.num))
py.readlast(verbose = False)
compareNumpyArraysClose(zeros, py.p_f, "Conservation of pressure:",

# Import sphere functionality
import sphere

initialization = True
plots = True

# CUDA device to use
device = 0

# Number of particles
np = 1e4

# Common simulation id
sim_id = 'init-1e4'

init = sphere.sim(np=np, nd=3, nw=0, sid=sim_id)

# Save radii
init.generateRadii(mean=0.01)

# Use default params
init.defaultParams(gamma_n=100.0, mu_s=0.6, mu_d=0.6)
init.setStiffnessNormal(1.16e7)
init.setStiffnessTangential(1.16e7)

# Add gravity
init.g[2] = -9.81

# Periodic x and y boundaries
init.periodicBoundariesXY()

#!/usr/bin/env python
from pytestutils import *
import sphere
import numpy

#'''
print("### Steady state, no gravity, no forcing, Dirichlet+Dirichlet BCs")
orig = sphere.sim('darcy_particles', np = 1000)
orig.cleanup()
#orig.generateRadii(histogram=False, psd='uni', radius_mean=5.0e-4, radius_variance=5.0e-5)
orig.defaultParams()
orig.generateRadii(psd='uni', mean=5.0e-2, variance=5.0e-5)
orig.initRandomGridPos([20, 20, 200])
orig.initTemporal(total=0.005, file_dt=0.001)
orig.initFluid(cfd_solver=1)
#orig.p_f[5,3,2] *= 1.5
#orig.k_c[0] = 4.6e-15
orig.k_c[0] = 4.6e-10
#orig.g[2] = -10.0
orig.setStiffnessNormal(36.4e9)
orig.setStiffnessTangential(36.4e9/3.0)
orig.run(verbose=False)
#orig.writeVTKall()
py = sphere.sim(sid = orig.sid, fluid = True)
py.readlast(verbose=False)

zeros = numpy.zeros((orig.num))
py.readlast(verbose = False)
compareNumpyArrays(zeros, py.p_f, "Conservation of pressure:")

# Fluid flow should be very small

#!/usr/bin/env python
from pytestutils import *

import sphere
#import sys
import numpy
#import matplotlib.pyplot as plt

print("### CFD tests - Dirichlet BCs ###")

# Iteration and conservation of mass test
# No gravity, no pressure gradients => no flow
print("# No forcing")
orig = sphere.sim(np = 0, nd = 3, nw = 0, sid = "cfdtest", fluid = True)
cleanup(orig)
orig.defaultParams()
#orig.defineWorldBoundaries([1.0,1.0,1.0], dx=0.1)
orig.defineWorldBoundaries([0.4,0.3,0.4], dx=0.1)
orig.initFluid(cfd_solver = 1)
#orig.initFluid(mu = 8.9e-4)
orig.initTemporal(total = 0.2, file_dt = 0.01, dt = 1.0e-7)
#orig.g[2] = -10.0
orig.time_file_dt = orig.time_dt*0.99
orig.time_total = orig.time_dt*10
#orig.run(dry=True)
py = sphere.sim(sid = orig.sid, fluid = True)
orig.run(verbose=False)
#orig.run(verbose=True)

zeros = numpy.zeros((orig.num))
py.readlast(verbose = False)

#!/usr/bin/env python
import sphere
from pytestutils import *

sim = sphere.sim('fluid_particle_interaction', fluid=True)
sim.cleanup()

sim.defineWorldBoundaries([1.0, 1.0, 1.0], dx = 0.1)
sim.initFluid()
sim.rho[0] = 1000.0 # particle density = fluid density
sim.setDEMstepsPerCFDstep(100)


# No gravity, pressure gradient enforced by Dirichlet boundaries.
# The particle should be sucked towards the low pressure
print('# Test 1: Test pressure gradient force')
sim.p_f[:,:,0]  = 10.0
sim.p_f[:,:,-1] = 1.0
sim.addParticle([0.5, 0.5, 0.5], 0.05)
sim.initTemporal(total=0.01, file_dt=0.001)

sim.run(verbose=False)
#sim.run(dry=True)
#sim.run(cudamemcheck=True)
#sim.writeVTKall()

sim.readlast()
test(sim.vel[0,2] > 0.0, 'Particle velocity:')

#!/usr/bin/env python
import sphere
from pytestutils import *

orig = sphere.sim('cfd_incl', fluid=True)
orig.cleanup()
#orig.defineWorldBoundaries([0.3, 0.3, 0.3], dx = 0.1)
orig.defineWorldBoundaries([0.3, 0.3, 0.3], dx = 0.06)
orig.initFluid(mu=8.9e-4) # inviscid "fluids" (mu=0) won't work!
#orig.initTemporal(total = 0.5, file_dt = 0.05, dt = 1.0e-4)
orig.initTemporal(total = 1.0e-0, file_dt = 1.0e-1, dt = 1.0e-3)
orig.bc_bot[0] = 1 # No-flow, free slip BC at bottom (Neumann)
#orig.bc_bot[0] = 2 # No-flow, no slip BC at bottom (Neumann)
#orig.bc_top[0] = 1 # No-flow, free slip BC at top (Neumann)

angle = 10.0 # slab inclination in degrees
g_magnitude = 10.0
orig.g[0] =  numpy.sin(numpy.radians(angle))*g_magnitude
orig.g[2] = -numpy.cos(numpy.radians(angle))*g_magnitude

tau_d = orig.g * orig.rho_f * orig.L[2] # analytical driving stress
v_sur = tau_d * orig.L[2] / orig.mu     # analytical surface velocity

# increase the max iterations for first step
orig.setMaxIterations(1e5)

# Homogeneous pressure, no gravity
orig.run(verbose=False)
orig.writeVTKall()

py = sphere.sim(sid=orig.sid, fluid=True)

device = int(sys.argv[1])
wet = int(sys.argv[2])
c_phi = float(sys.argv[3])
c_grad_p = float(sys.argv[4])
sigma0 = float(sys.argv[5])

#sim = sphere.sim('diffusivity-sigma0=' + str(sigma0) + '-c_phi=' + \
#        str(c_phi) + '-c_grad_p=' + str(c_grad_p), fluid=True)
if wet == 1:
    fluid = True
else:
    fluid = False
    
#sim = sphere.sim('diffusivity-sigma0=' + str(sigma0) +'-c_phi=1.0-c_grad_p=1.0',
#        fluid=True)
sim = sphere.sim('diffusivity-sigma0=' + str(sigma0),
        fluid=True)
sim.readlast()

if sigma0 == 20.0e3 and c_phi == 1.0 and c_grad_p == 0.1:
    sim.sid = 'shear-sigma0=20000.0-c_phi=1.0-c_grad_p=0.1-hi_mu-lo_visc-hw-noshear'
    sim.readlast()

if fluid:
    sim.sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + str(c_phi) + \
            '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc-hw'
else:
    sim.sid = 'shear-sigma0=' + str(sigma0) + '-hw'

print(sim.sid)
sim.fluid = fluid