Skip to content

Instantly share code, notes, and snippets.

@wd15

wd15/README.md

Last active May 17, 2017 14:46
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save wd15/c28ab796cb3d9781482b01fb67a7ec2d to your computer and use it in GitHub Desktop.
Save wd15/c28ab796cb3d9781482b01fb67a7ec2d to your computer and use it in GitHub Desktop.
Download ZIP
Reactive Wetting Code
Raw
README.md

Reactive Wetting Using FiPy

Raw
density.py
from fipy import *
from dictTable import DictTable
from triplePoint import TriplePointContour
from matplotlib import rc, rcParams
import matplotlib
matplotlib.use('Agg')
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rc('text', usetex=True)
import pylab
pylab.ioff()
pylab.hold(True)
Ly = 3.6
Lx = 2.5
size = 4.
pylab.figure(figsize=(2 * Lx / Ly * size, size))
dataFile = os.path.join('no_tsm', 'data.h5')
data = DictTable(dataFile, 'r')
print data.getLatestIndex()
for timestep, color in ((0, 'b'),
## (974, 'g'),
## (984, 'g'),
(data.getLatestIndex(), 'r')):
## (2000, 'r'),
## (2200, 'r'),
## (2400, 'r'),
## (2600, 'r'),
## (2800, 'r'),
## (3000, 'r'),
## (3200, 'r')):
dx = data[timestep]['dx']
dy = data[timestep]['dy']
nx = data[timestep]['nx']
ny = data[timestep]['ny']
mesh = CylindricalGrid2D(nx=nx, dx=dx, ny=ny, dy=dy)
density1 = CellVariable(mesh=mesh, value=data[timestep]['density1'])
density2 = CellVariable(mesh=mesh, value=data[timestep]['density2'])
Xvelocity = CellVariable(mesh=mesh, value=data[timestep]['Yvelocity'])
Yvelocity = CellVariable(mesh=mesh, value=data[timestep]['Xvelocity'])
density = density1 + density2
phase = CellVariable(mesh=mesh, value=data[timestep]['phase'])
X, Y = mesh.getCellCenters()
print
print 'X[0]',X[0]
print 'mesh.getCellVolumes()[0]',mesh.getCellVolumes()[0]
print 'mesh.getCellVolumes()[1]',mesh.getCellVolumes()[1]
print
print 'density([(X[0],), (2.5,)])',density([(X[0],), (2.5e-6,)])
print 'density1([(X[0],), (2.5,)])',density1([(X[0],), (2.5e-6,)])
print 'density2([(X[0],), (2.5,)])',density2([(X[0],), (2.5e-6,)])
print 'Xvelocity([(X[0],), (2.5,)])',Xvelocity([(X[0],), (2.5,)])
print 'Yvelocity([(X[0],), (2.5,)])',Yvelocity([(X[0],), (2.5,)])
print
print 'density([(X[1],), (2.5,)])',density([(X[1],), (2.5e-6,)])
print 'density1([(X[1],), (2.5,)])',density1([(X[1],), (2.5e-6,)])
print 'density2([(X[1],), (2.5,)])',density2([(X[1],), (2.5e-6,)])
print 'Xvelocity([(X[1],), (2.5,)])',Xvelocity([(X[1],), (2.5,)])
print 'Yvelocity([(X[1],), (2.5,)])',Yvelocity([(X[1],), (2.5,)])
ld1=767.671770874
ld2=6586.6684953599997
sd1=7489.1017183100003
sd2=349.289285769
vd1=8.64877867561
vd2=74.207024652200005
vd = vd1 + vd2
ld = ld1 + ld2
phase = (phase > 0.) * phase + (phase > 1.) * ( 1. - phase)
gamma = (density - vd) / (ld - vd)
gamma = (gamma > 0.) * gamma + (gamma > 1.) * ( 1. - gamma)
phis = phase
phil = gamma * (1 - phase)
phiv = (1 - gamma) * (1 - phase)
##nx = N
##ny = N
##nx = int((Xrange[1] - Xrange[0]) / dx + 1e-10)
##ny = int((Yrange[1] - Yrange[0]) / dy + 1e-10)
phis = numerix.reshape(phis, (ny, nx))
phil = numerix.reshape(phil, (ny, nx))
tmps = numerix.zeros((ny, 2 * nx), 'd')
tmpl = numerix.zeros((ny, 2 * nx), 'd')
tmps[:,:nx] = phis[:,::-1]
tmps[:,nx:] = phis
tmpl[:,:nx] = phil[:,::-1]
tmpl[:,nx:] = phil
phis = tmps
phil = tmpl
pylab.contour(phis, (0.5,), extent = (-dx * nx * 1e+6, dx * nx * 1e+6, -0.5, dy * ny * 1e+6 - 0.5), colors=color)
pylab.contour(phil, (0.5,), extent = (-dx * nx * 1e+6, dx * nx * 1e+6, -0.5, dy * ny * 1e+6 - 0.5), colors=color)
R0 = 1.0
solidHeight = 0.5
h0 = Ly / 2. - solidHeight
gamma_lv = 18.98
gamma_sv = 61.14 - 18.98
gamma_sl = 31.74
theta = numerix.arccos((gamma_sv - gamma_sl) / gamma_lv)
def Vcutt(R, h):
hbar = R - h
return numerix.pi * R * hbar**2 - numerix.pi * hbar**3 / 3.
initialVolume = 4. * numerix.pi * R0**3 / 3. ##- Vcutt(R0, h0)
div = sqrt(1. + tan(theta)**2)
R1 = (initialVolume / Vcutt(1., 1. / div))**(1. / 3.)
h1 = sqrt(R1**2 / div**2)
r1 = h1 * tan(theta)
pylab.plot((-Ly, Ly), (0, 0), 'k')
N = 200
x = 2. * (arange(N) * r1 / (N - 1)) - r1
y = numerix.sqrt(R1**2 - x**2) - h1
print 'finalVolume',Vcutt(R1, h1)
pylab.plot(x , y, 'k')
pylab.xlim(xmin=-2.5, xmax=2.5)
pylab.xlabel(r'$r$ ($\mu$m)')
pylab.ylabel(r'$h$ ($\mu$m)')
pylab.ylim(ymin=-solidHeight, ymax=Ly - solidHeight)
pylab.savefig('density.pdf')
Raw
density1.py
from fipy import *
from dictTable import DictTable
from triplePoint import TriplePointContour
from matplotlib import rc, rcParams
import matplotlib
matplotlib.use('Agg')
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rc('text', usetex=True)
import pylab
pylab.ioff()
pylab.hold(True)
Ly = 3.6
Lx = 2.5
size = 4.
pylab.figure(figsize=(2 * Lx / Ly * size, size))
dataFile = os.path.join('no_tsm', 'data.h5')
data = DictTable(dataFile, 'r')
print data.getLatestIndex()
for timestep, color in ((0, 'k'),
(500, 'k'),
(1000, 'k'),
(1500, 'k'),
(2000, 'k'),
(2500, 'k'),
(3000, 'k'),
(data.getLatestIndex(), 'k')):
## (2000, 'r'),
## (2200, 'r'),
## (2400, 'r'),
## (2600, 'r'),
## (2800, 'r'),
## (3000, 'r'),
## (3200, 'r')):
dx = data[timestep]['dx']
dy = data[timestep]['dy']
nx = data[timestep]['nx']
ny = data[timestep]['ny']
mesh = CylindricalGrid2D(nx=nx, dx=dx, ny=ny, dy=dy)
density1 = CellVariable(mesh=mesh, value=data[timestep]['density1'])
density2 = CellVariable(mesh=mesh, value=data[timestep]['density2'])
density = density1 + density2
phase = CellVariable(mesh=mesh, value=data[timestep]['phase'])
X, Y = mesh.getCellCenters()
ld1=767.671770874
ld2=6586.6684953599997
sd1=7489.1017183100003
sd2=349.289285769
vd1=8.64877867561
vd2=74.207024652200005
vd = vd1 + vd2
ld = ld1 + ld2
phase = (phase > 0.) * phase + (phase > 1.) * ( 1. - phase)
gamma = (density - vd) / (ld - vd)
gamma = (gamma > 0.) * gamma + (gamma > 1.) * ( 1. - gamma)
phis = phase
phil = gamma * (1 - phase)
phiv = (1 - gamma) * (1 - phase)
##nx = N
##ny = N
##nx = int((Xrange[1] - Xrange[0]) / dx + 1e-10)
##ny = int((Yrange[1] - Yrange[0]) / dy + 1e-10)
phis = numerix.reshape(phis, (ny, nx))
phil = numerix.reshape(phil, (ny, nx))
tmps = numerix.zeros((ny, 2 * nx), 'd')
tmpl = numerix.zeros((ny, 2 * nx), 'd')
tmps[:,:nx] = phis[:,::-1]
tmps[:,nx:] = phis
tmpl[:,:nx] = phil[:,::-1]
tmpl[:,nx:] = phil
phis = tmps
phil = tmpl
pylab.contour(phis, (0.5,), extent = (-dx * nx * 1e+6, dx * nx * 1e+6, -0.5, dy * ny * 1e+6 - 0.5), colors=color, linewidths=0.4)
pylab.contour(phil, (0.5,), extent = (-dx * nx * 1e+6, dx * nx * 1e+6, -0.5, dy * ny * 1e+6 - 0.5), colors=color, linewidths=0.4)
R0 = 1.0
solidHeight = 0.5
h0 = Ly / 2. - solidHeight
gamma_lv = 18.98
gamma_sv = 61.14 - 18.98
gamma_sl = 31.74
theta = numerix.arccos((gamma_sv - gamma_sl) / gamma_lv)
def Vcutt(R, h):
hbar = R - h
return numerix.pi * R * hbar**2 - numerix.pi * hbar**3 / 3.
initialVolume = 4. * numerix.pi * R0**3 / 3. ##- Vcutt(R0, h0)
div = sqrt(1. + tan(theta)**2)
R1 = (initialVolume / Vcutt(1., 1. / div))**(1. / 3.)
h1 = sqrt(R1**2 / div**2)
r1 = h1 * tan(theta)
##pylab.plot((-Ly, Ly), (0, 0), 'k')
N = 200
x = 2. * (arange(N) * r1 / (N - 1)) - r1
y = numerix.sqrt(R1**2 - x**2) - h1
print 'finalVolume',Vcutt(R1, h1)
##pylab.plot(x , y, 'k')
pylab.xlim(xmin=-2.5, xmax=2.5)
pylab.xlabel(r'$r$ ($\mu$m)')
pylab.ylabel(r'$h$ ($\mu$m)')
pylab.ylim(ymin=-solidHeight, ymax=Ly - solidHeight)
pylab.savefig('density1.pdf')
Raw
dictTable.py
from fipy import *
import tables
import os
class DictTable:
"""
Designed to save a dictionary of arrays at each timestep in a simulation.
"""
IDprefix = 'ID'
def __init__(self, h5filename, mode='a'):
if mode == 'w' and os.path.exists(h5filename):
print 'removing %s ' % h5filename
os.remove(h5filename)
self.h5filename = h5filename
## self.h5file = tables.openFile(h5filename, mode=mode)
def __setitem__(self, index, values):
h5file = tables.openFile(self.h5filename, mode='a')
## self.h5file.mode = 'w'
## h5file.mode = 'w'
## group = self.h5file.createGroup(self.h5file.root, self.IDprefix + str(index))
h5file.root._v_attrs.latestIndex = index
groupName = self.IDprefix + str(index)
if hasattr(h5file.root, groupName):
group = h5file.root._f_getChild(groupName)
group._f_remove(recursive=True)
group = h5file.createGroup(h5file.root, groupName)
for k in values.keys():
## self.h5file.createArray(group, k, values[k])
h5file.createArray(group, k, values[k])
h5file.close()
def __getitem__(self, index):
## self.h5file.mode = 'r'
h5file = tables.openFile(self.h5filename, mode='r')
d = {}
## for array in self.h5file.listNodes('/' + self.IDprefix + str(index), classname='Array'):
for array in h5file.listNodes('/' + self.IDprefix + str(index), classname='Array'):
d[array.name] = array.read()
h5file.close()
return d
def getLatestIndex(self):
h5file = tables.openFile(self.h5filename, mode='r')
latestIndex = h5file.root._v_attrs.latestIndex
h5file.close()
return latestIndex
Raw
dropletShape.py
## Find file will evaluate the radius of a drop.
## Assumptions: flat solid vapor interface, drop equilibrium densities are not adjusted
def dropletShape(gamma_lv=None, gamma_sv=None, gamma_sl=None, L=None, ld1=None, ld2=None, sd1=None, sd2=None, vd1=None, vd2=None):
## evaluate the upper and lower angles
c_upper = (gamma_sv**2 + gamma_lv**2 - gamma_sl**2) / (2 * gamma_sv * gamma_lv)
c_lower = (gamma_sv**2 - gamma_lv**2 + gamma_sl**2) / (2 * gamma_sv * gamma_sl)
import numpy
s_upper = numpy.sqrt(1 - c_upper**2)
s_lower = numpy.sqrt(1 - c_lower**2)
theta_upper = numpy.arccos(c_upper)
theta_lower = numpy.arccos(c_lower)
## print 'theta_upper',theta_upper * 180. / numpy.pi
## print 'theta_lower',theta_lower * 180. / numpy.pi
## check to see that the values for the cosines satisfy the law of sines.
## print 'gamma_lv * c_upper + gamma_sl * c_lower - gamma_sv = 0 =',gamma_lv * c_upper + gamma_sl * c_lower - gamma_sv
## print 'gamma_sl * s_lower - gamma_lv * s_upper = 0 =',gamma_sl * s_lower - gamma_lv * s_upper
## the initial configuration
initialSolidHeight = L / 4.
## widthOfBox = 3. * L / 4.
widthOfBox = L
initialRadius = L / 2.
vd1_eq = vd1
ld1_eq = ld1
sd1_eq = sd1
vd2_eq = vd2
ld2_eq = ld2
sd2_eq = sd2
vd1_initial = vd1_eq
ld1_initial = ld1_eq
sd1_initial = sd1_eq
vd2_initial = vd2_eq
ld2_initial = ld2_eq
sd2_initial = sd2_eq
Q1 = widthOfBox * initialSolidHeight * (sd1_initial - vd1_initial) \
+ numpy.pi * initialRadius**2 / 4. * (ld1_initial - vd1_initial) + widthOfBox * L * vd1_initial
Q2 = widthOfBox * initialSolidHeight * (sd2_initial - vd2_initial) \
+ numpy.pi * initialRadius**2 / 4. * (ld2_initial - vd2_initial) + widthOfBox * L * vd2_initial
## final configuration
vd1_final = vd1_eq
ld1_final = ld1_eq
sd1_final = sd1_eq
vd2_final = vd2_eq
ld2_final = ld2_eq
sd2_final = sd2_eq
A_bar_upper = (theta_upper - c_upper * s_upper) / s_upper**2 / 2.
A_bar_lower = (theta_lower - c_lower * s_lower) / s_lower**2 / 2.
A_bar = A_bar_upper + A_bar_lower
a1 = A_bar * ld1_final - A_bar_lower * sd1_final - A_bar_upper * vd1_final
a2 = A_bar * ld2_final - A_bar_lower * sd2_final - A_bar_upper * vd2_final
b1 = (sd1_final - vd1_final) * widthOfBox
b2 = (sd2_final - vd2_final) * widthOfBox
h = (Q1 * a2 - Q2 * a1) / (b1 * a2 - a1 * b2)
R = numpy.sqrt((Q1 * b2 - Q2 * b1) / (a1 * b2 - b1 * a2))
return h, R, theta_upper, theta_lower
print 'final height of solid is',h
print 'radius of drop is',R
if __name__ == '__main__':
gamma_lv = (18.148 + 18.852) / 2.
gamma_sv = (24.572 + 20.135) / 2.
gamma_sl = (26.014 + 25.931) / 2.
L = 1e-6
ld1 = 2944.01311655
sd1 = 6215.54674682
ld2 = 4325.00238504
sd2 = 1867.77409607
h, R, theta1, theta2 = dropletShape(gamma_lv=gamma_lv,
gamma_sv=gamma_sv,
gamma_sl=gamma_sl,
L=L,
ld1=ld1,
ld2=ld2,
sd1=sd1,
sd2=sd2)
print 'final height of solid is',h
print 'radius of drop is',R
Raw
input.py
#!/usr/bin/env python
from fipy import *
factor = 1
Lx = 2.5e-6 ##* factor
Ly = 3.6e-6 ##* factor
nx = 250 * factor
ny = 360 * factor
from PyTrilinos import Epetra
Nproc = Epetra.PyComm().NumProc()
print Epetra.PyComm().MyPID()
print Nproc
ny = ny + (Nproc - ny % Nproc) * ((ny % Nproc) > 0)
dx = Lx / nx
dy = Ly / ny
viscosity_liquid = 2e-2
viscosity_solid = 2e+6
viscosity_vapor = 2e-2
epsilon1 = 1e-16
epsilon2 = 1e-16
temperature = 650.0
molarWeight = 0.118
gasConstant = 8.314
totalTimeSteps = 100000000000000
cfl = 0.1
va = 1.0
tolerance = 1e-1
ee = -0.455971
vbar = 1.3e-05
initial_dt = 1e-15
dt = initial_dt
A1 = 28263.829467072799
A2 = 56443.954187683703
B = 25418.658283299999 * 4
vs = 1.50702137711749e-05
W = 127408.43535226236
epsilon_phi = 1.0e-09
mobility_phi = 10000.0
density1Solid = 7489.1017183100003
density2Solid = 349.289285769
density1Liquid = 767.671770874
density1Vapor = 8.64877867561
density2Liquid = 6586.6684953599997
density2Vapor = 74.207024652200005
densitySolid = density1Solid + density2Solid
densityLiquid = density1Liquid + density2Liquid
densityVapor = density1Vapor + density2Vapor
delta_conc = 1.0
density1LiquidNew = delta_conc * density1Liquid
density2LiquidNew = densityLiquid - delta_conc * density1Liquid
density1Liquid = density1LiquidNew
density2Liquid = density2LiquidNew
max_sweeps = 20
Mbar_fluid = 1e-7
Mbar_solid= 1e-11
##from fipy.meshes.cylindricalGrid2D import CylindricalGrid2D
from parallelCylindricalGrid2D import ParallelCylindricalGrid2D
shift = dx * 1e-10
mesh = ParallelCylindricalGrid2D(nx=nx, ny=ny, dx=dx, dy=dy, origin=((shift,),(0,)))
density1 = CellVariable(mesh=mesh, hasOld=True)
density2 = CellVariable(mesh=mesh, hasOld=True)
phaseField = CellVariable(mesh=mesh, hasOld=True)
## initial position ##
X, Y = mesh.getCellCenters()
X = X - shift
Xf, Yf = mesh.getFaceCenters()
Xf = Xf - shift
phaseField.setValue(0)
density1.setValue(density1Vapor)
density2.setValue(density2Vapor)
solidHeight = 0.5e-6
dropRadius = 1e-6
dropCenter = (shift, dropRadius + solidHeight)
dropMask = (X - dropCenter[0])**2 + (Y - dropCenter[1])**2 < dropRadius**2
phaseField.setValue(0, where=dropMask)
density1.setValue(density1Liquid, where=dropMask)
density2.setValue(density2Liquid, where=dropMask)
phaseField.setValue(1, where=Y < (solidHeight - 1e-20))
density1.setValue(density1Solid, where=Y < (solidHeight - 1e-20))
density2.setValue(density2Solid, where=Y < (solidHeight - 1e-20))
phaseBoundary = 1.0
density1.updateOld()
density2.updateOld()
phaseField.updateOld()
####################### smoothing #####################
from unsegregatedEquation import UnsegregatedEquation
smoothPhaseEqn = UnsegregatedEquation(((TransientTerm(1) - DiffusionTerm(coeff=1),),))
smoothDensityEqn = UnsegregatedEquation(((TransientTerm(1) - DiffusionTerm(coeff=1),),))
smoothX1Eqn = UnsegregatedEquation(((TransientTerm(1) - DiffusionTerm(coeff=1),),))
fakedt = 1e-15
smoothDensity = CellVariable(mesh=mesh, hasOld=False)
smoothX1 = CellVariable(mesh=mesh, hasOld=False)
smoothDensity.setValue(density1 + density2)
smoothX1.setValue(density1 / (density1 + density2))
smoothPhaseEqn.sweep((phaseField,), dt=fakedt, boundaryConditions=(((),),), relaxation=1.0)
smoothDensityEqn.sweep((smoothDensity,), dt=fakedt, boundaryConditions=(((),),), relaxation=1.0)
smoothDensityEqn.sweep((smoothX1,), dt=fakedt, boundaryConditions=(((),),), relaxation=1.0)
density1.setValue(smoothX1 * smoothDensity)
density2.setValue((1 - smoothX1) * smoothDensity)
del smoothDensity
del smoothX1
########################################################
Xvelocity = CellVariable(mesh=mesh, hasOld=True)
Yvelocity = CellVariable(mesh=mesh, hasOld=True)
NC1 = CellVariable(mesh=mesh, hasOld=False)
NC2 = CellVariable(mesh=mesh, hasOld=False)
faceVelocity = FaceVariable(mesh=mesh, rank=1)
ap = CellVariable(mesh=mesh, value=1e+20)
density = density1 + density2
faceDensity = density.getArithmeticFaceValue()
faceDensity1 = density1.getArithmeticFaceValue()
faceDensity2 = density2.getArithmeticFaceValue()
X1_face = faceDensity1 / faceDensity
X2_face = faceDensity2 / faceDensity
phaseFieldSq = phaseField**2
phaseFieldSqSq = phaseFieldSq**2
Mbar = Mbar_solid**phaseFieldSqSq * Mbar_fluid**(1 - phaseFieldSqSq)
M = Mbar.getHarmonicFaceValue() * X1_face * X2_face
alpha0 = (density - densityVapor) / (densityLiquid - densityVapor)
alpha = (alpha0 > 0.) * alpha0 + (alpha0 > 1.) * ( 1. - alpha0)
viscosity_fluid = viscosity_liquid**alpha * viscosity_vapor**(1 - alpha)
viscosity = (viscosity_solid**phaseFieldSqSq * viscosity_fluid**(1 - phaseFieldSqSq)).getMinmodFaceValue()
ConvectionTerm = CentralDifferenceConvectionTerm
## liquid chemical potentials
div = molarWeight - vbar * density
RTM = gasConstant * temperature / molarWeight
RTMdivsq = RTM / div**2
twoeeM = 2 * ee / molarWeight**2
pp1 = vbar * (vbar * density - 2 * molarWeight)
pp2 = molarWeight * vbar * density / div**2
mu1_derivative1_l = twoeeM + RTMdivsq / density1 * (molarWeight**2 + density2 * pp1)
mu1_derivative2_l = twoeeM - RTMdivsq * pp1
mu1_const_l = RTM * (numerix.log(density1 / div) - pp2)
mu2_derivative2_l = twoeeM + RTMdivsq / density2 * (molarWeight**2 + density1 * pp1)
mu2_derivative1_l = mu1_derivative2_l
mu2_const_l = RTM * (numerix.log(density2 / div) - pp2)
## solid chemcial potentials
pp4 = B * molarWeight / density**3 / vs**2
pp3 = -RTM / density + 2 * pp4
pp5 = B / molarWeight - 3 * B * molarWeight / density**2 / vs**2
mu1_derivative1_s = RTM / density1 + pp3
mu1_derivative2_s = pp3
mu1_const_s = A1 / molarWeight + RTM * log(density1 / density) + pp5
mu2_derivative2_s = RTM / density2 + pp3
mu2_derivative1_s = mu1_derivative2_s
mu2_const_s = A2 / molarWeight + RTM * log(density2 / density) + pp5
## chemical potentials
p = phaseField**3 * (10 - 15 * phaseField + 6 * phaseField**2)
pPrime = 30 * phaseField**2 * (1 - phaseField)**2
pPrimePrime = 60 * phaseField * (1 - phaseField) * (1 - 2 * phaseField)
oneMinusP = 1 - p
mu1_derivative1 = p * mu1_derivative1_s + oneMinusP * mu1_derivative1_l
mu1_derivative2 = p * mu1_derivative2_s + oneMinusP * mu1_derivative2_l
mu1_const = p * mu1_const_s + oneMinusP * mu1_const_l
mu2_derivative2 = p * mu2_derivative2_s + oneMinusP * mu2_derivative2_l
mu2_derivative1 = p * mu2_derivative1_s + oneMinusP * mu2_derivative1_l
mu2_const = p * mu2_const_s + oneMinusP * mu2_const_l
## phase potential
fsminusfl = A1 * density1 / molarWeight + A2 * density2 / molarWeight \
+ B * molarWeight / density / vs**2 * (1 - density * vs / molarWeight)**2 \
- RTM * density * numerix.log(density / div) \
- ee * density**2 / molarWeight**2
dfdphi = fsminusfl * pPrime + W / 30 * pPrimePrime
## velocity correction
apCoeff = mesh._getFaceAreas() * mesh._getCellDistances() / ap.getArithmeticFaceValue()
density1Coeff = apCoeff * faceDensity1
density2Coeff = apCoeff * faceDensity2
NCphasePotential = dfdphi - epsilon_phi* temperature * phaseField.getFaceGrad().getDivergence()
tensorDiff = - phaseField.getFaceGrad() * NCphasePotential.getArithmeticFaceValue()
tensorAverage = (density1 * NC1.getLeastSquaresGrad() + density2 * NC2.getLeastSquaresGrad() - phaseField.getLeastSquaresGrad() * NCphasePotential).getArithmeticFaceValue()
velocityCorrection = apCoeff * (tensorDiff - tensorAverage)
from unsegregatedEquation import UnsegregatedEquation
Xface, Yface = mesh.getFaceCenters()
uu = FaceVariable(mesh=mesh, rank=1)
uu[0] = 1 / Xface
normal_x = (1, 0)
matx = ((0, 0), (1, 0))
normal_y = (0, 1)
maty = ((0, 1), (0, 0))
combinedEquation = UnsegregatedEquation(((TransientTerm(density) \
+ ConvectionTerm(faceDensity * faceVelocity) \
- DiffusionTerm((viscosity + normal_x * viscosity,)),
- DiffusionTerm((matx * viscosity,)),
None,
None,
X * density1 * ConvectionTerm(uu),
X * density2 * ConvectionTerm(uu),
phaseField.getLeastSquaresGrad()[0] * (dfdphi - DiffusionTerm(epsilon_phi * temperature))),
(- DiffusionTerm((maty * viscosity,)),
TransientTerm(density) \
+ ConvectionTerm(faceDensity * faceVelocity) \
- DiffusionTerm((viscosity + normal_y * viscosity,)),
None,
None,
density1 * ConvectionTerm((0, 1)),
density2 * ConvectionTerm((0, 1)),
phaseField.getLeastSquaresGrad()[1] * (dfdphi - DiffusionTerm(epsilon_phi * temperature))),
(None,
None,
TransientTerm() \
+ VanLeerConvectionTerm(faceVelocity - velocityCorrection),
None,
DiffusionTerm(-M / temperature - density1Coeff * faceDensity1),
+ DiffusionTerm(M / temperature - density1Coeff * faceDensity2),
None),
(None,
None,
None,
TransientTerm() \
+ VanLeerConvectionTerm(faceVelocity - velocityCorrection),
DiffusionTerm(M / temperature - density2Coeff * faceDensity1),
+ DiffusionTerm(-M / temperature - density2Coeff * faceDensity2),
None),
(None,
None,
-ImplicitSourceTerm(mu1_derivative1) + DiffusionTerm(epsilon1 * temperature),
-ImplicitSourceTerm(mu1_derivative2),
ImplicitSourceTerm(1.) - mu1_const,
None,
None),
(None,
None,
-ImplicitSourceTerm(mu2_derivative1),
-ImplicitSourceTerm(mu2_derivative2) + DiffusionTerm(epsilon2 * temperature),
None,
ImplicitSourceTerm(1.) - mu2_const,
None),
(None,
None,
None,
None,
None,
None,
TransientTerm() + X * Xvelocity * ConvectionTerm(uu) + Yvelocity * ConvectionTerm((0,1)) \
- DiffusionTerm(epsilon_phi * mobility_phi) \
+ mobility_phi / temperature * dfdphi)))
bc00 = (FixedValue(value=0., faces=mesh.getFacesLeft() + mesh.getFacesRight()),
FixedFlux(value=0., faces=mesh.getFacesUp() + mesh.getFacesDown()))
bc04 = (FixedValue(value=0., faces=mesh.getExteriorFaces()), FixedValue(value=0., faces=mesh.getExteriorFaces()))
bc05 = (FixedValue(value=0., faces=mesh.getExteriorFaces()), FixedValue(value=0., faces=mesh.getExteriorFaces()))
bc11 = (FixedFlux(value=0., faces=mesh.getFacesLeft() + mesh.getFacesRight()),
FixedValue(value=0., faces=mesh.getFacesUp() + mesh.getFacesDown()))
bc14 = (FixedValue(value=0., faces=mesh.getExteriorFaces()), FixedValue(value=0., faces=mesh.getExteriorFaces()))
bc15 = (FixedValue(value=0., faces=mesh.getExteriorFaces()), FixedValue(value=0., faces=mesh.getExteriorFaces()))
bc66 = (FixedValue(value=phaseBoundary, faces=mesh.getFacesBottom()),)
## exteriorFaces = mesh.getFacesUp() + mesh.getFacesDown() + mesh.getFacesRight()
## bc00 = (FixedValue(value=0., faces=mesh.getFacesRight()),
## FixedFlux(value=0., faces=mesh.getFacesUp() + mesh.getFacesDown()))
## bc04 = (FixedValue(value=0., faces=exteriorFaces), FixedValue(value=0., faces=exteriorFaces))
## bc05 = (FixedValue(value=0., faces=exteriorFaces), FixedValue(value=0., faces=exteriorFaces))
## bc11 = (FixedFlux(value=0., faces=mesh.getFacesRight()),
## FixedValue(value=0., faces=mesh.getFacesUp() + mesh.getFacesDown()))
## bc14 = (FixedValue(value=0., faces=mesh.getExteriorFaces()), FixedValue(value=0., faces=mesh.getExteriorFaces()))
## bc15 = (FixedValue(value=0., faces=mesh.getExteriorFaces()), FixedValue(value=0., faces=mesh.getExteriorFaces()))
## bc66 = (FixedValue(value=phaseBoundary, faces=mesh.getFacesBottom()),)
BCs = ((bc00, () , (), (), bc04, bc05, () ),
(() , bc11, (), (), bc14, bc15, () ),
(() , () , (), (), () , () , () ),
(() , () , (), (), () , () , () ),
(() , () , (), (), () , () , () ),
(() , () , (), (), () , () , () ),
(() , () , (), (), () , () , bc66))
from dictTable import DictTable
from PyTrilinos import Epetra
##procID = Epetra.PyComm().MyPID()
import os.path
dataFile = os.path.join('no_tsm', 'data.h5')
def aggregateData(numpyVector):
IDs = mesh.getGlobalCellIDs()
comm = Epetra.PyComm()
localMap = Epetra.Map(-1, list(IDs), 0, comm)
localVector = Epetra.Vector(localMap, numpyVector)
globalMap = Epetra.Map(-1, list(numerix.arange(nx * ny)), 0, comm)
globalVector = Epetra.Vector(globalMap)
globalVector.Import(localVector, Epetra.Import(globalMap, localMap), Epetra.Insert)
return numerix.array(globalVector)
def writeData(h5data, timeStep):
dataDict = {'elapsedTime' : elapsedTime,
'timeStep' : int(timeStep),
'totalSweeps' : int(totalSweeps),
'dt' : dt,
'nx' : int(nx),
'ny' : int(ny),
'dx' : dx,
'dy' : dy,
'Nproc' : int(Epetra.PyComm().NumProc())}
for key, data in [('Xvelocity', numerix.array(Xvelocity)),
('Yvelocity', numerix.array(Yvelocity)),
('density1', numerix.array(density1)),
('density2', numerix.array(density2)),
('phase', numerix.array(phaseField)),
('ap', numerix.array(ap))]:
data = aggregateData(data)
if Epetra.PyComm().MyPID() == 0:
dataDict[key] = data
if Epetra.PyComm().MyPID() == 0:
h5data[int(timeStep)] = dataDict
if os.path.exists(dataFile):
from PyTrilinos import Epetra
procID = Epetra.PyComm().MyPID()
import os
for PID in range(Nproc):
if PID == procID:
print PID
## tmpfile = os.tempnam()
## import shutil
## shutil.copyfile(dataFile, tmpfile)
h5data = DictTable(dataFile, 'r')
## tmpdata = DictTable(tmpfile, 'r')
tmpdata = h5data
timeStep = tmpdata.getLatestIndex()
dt = tmpdata[timeStep]['dt']
totalSweeps = tmpdata[timeStep]['totalSweeps']
timeStep = tmpdata[timeStep]['timeStep']
elapsedTime = tmpdata[timeStep]['elapsedTime']
oldMesh = Grid2D(nx=tmpdata[timeStep]['nx'],
ny=tmpdata[timeStep]['ny'],
dx=tmpdata[timeStep]['dx'],
dy=tmpdata[timeStep]['dy'])
X, Y = mesh.getCellCenters()
for key, var in [('Xvelocity', Xvelocity),
('Yvelocity', Yvelocity),
('density1', density1),
('density2', density2),
('phase', phaseField),
('ap', ap)]:
oldVar = CellVariable(mesh=oldMesh)
oldVar.setValue(tmpdata[timeStep][key])
var.setValue(oldVar((X, Y), order=1))
## os.remove(tmpfile)
Epetra.PyComm().Barrier()
timeStep += 1
else:
timeStep = 0
totalSweeps = 0
elapsedTime = 0.0
if Epetra.PyComm().MyPID() == 0:
h5data = DictTable(dataFile, 'w')
else:
h5data = None
writeData(h5data, timeStep)
import time
tWall = time.time()
while True:
dt = min(dt * 1.1, dx / max(sqrt(Xvelocity**2 + Yvelocity**2)) * cfl)
dt = max(dt, initial_dt)
sweep = 0
residual = 1.
density1.updateOld()
density2.updateOld()
Xvelocity.updateOld()
Yvelocity.updateOld()
phaseField.updateOld()
NC1[:] = 0.0
NC2[:] = 0.0
while residual > tolerance:
dt = min(dt, dx / max(sqrt(Xvelocity**2 + Yvelocity**2)) * cfl)
dt = Epetra.PyComm().MinAll(dt)
faceVelocity[0] = Xvelocity.getArithmeticFaceValue()
faceVelocity[1] = Yvelocity.getArithmeticFaceValue()
faceVelocity[0, (Xf == 0) | (Xf == Lx)] = 0.
faceVelocity[1, (Yf == 0) | (Yf == Ly)] = 0.
previousResidual = residual
residual, matrixDiag, convergence = combinedEquation.sweep((Xvelocity,
Yvelocity,
density1,
density2,
NC1,
NC2,
phaseField), dt=dt, boundaryConditions=BCs)
if sweep < 2:
initialResidual = residual
residual = residual / initialResidual
redo_timestep = False
if residual > previousResidual or (convergence != 0 and convergence != -3) or sweep == max_sweeps:
redo_timestep = True
break
ap.setValue(matrixDiag[:mesh.getNumberOfCells()])
if Epetra.PyComm().MyPID() == 0:
print 'timeStep',timeStep,'dt',dt,'sweep',sweep,'residual',residual,'previousResidual',previousResidual,'convergence',convergence,'elapsedTime',elapsedTime,'wall time',time.time() - tWall
tWall = time.time()
totalSweeps += 1
sweep += 1
if redo_timestep is True:
if Epetra.PyComm().MyPID() == 0:
print
print 'residual',residual
print 'previousResidual',previousResidual
print 'convergence',convergence
print 'sweep',sweep
print 'max_sweeps',max_sweeps
density1.setValue(density1.getOld())
density2.setValue(density2.getOld())
Xvelocity.setValue(Xvelocity.getOld())
Yvelocity.setValue(Yvelocity.getOld())
phaseField.setValue(phaseField.getOld())
dt = dt / 5.
dt = max(dt , initial_dt)
else:
elapsedTime += dt
timeStep += 1
if timeStep % 10 == 0:
writeData(h5data, timeStep)
Raw
machinefile
benson.nist.gov:4
hobson.nist.gov:4
alice.nist.gov:4
Raw
movie.py
from fipy import *
from dictTable import DictTable
from triplePoint import TriplePointContour
from matplotlib import rc, rcParams
import matplotlib
import matplotlib.cm
print matplotlib.cm
matplotlib.use('Agg')
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rc('text', usetex=True)
import pylab
pylab.ioff()
pylab.hold(True)
Ly = 3.5
def makePNG(data, filename):
pylab.figure(figsize=(8.,4.))
dx = data['dx']
dy = data['dy']
nx = data['nx']
ny = data['ny']
mesh = CylindricalGrid2D(nx=nx, dx=dx, ny=ny, dy=dy)
density1 = CellVariable(mesh=mesh, value=data['density1'])
density2 = CellVariable(mesh=mesh, value=data['density2'])
density = density1 + density2
phase = CellVariable(mesh=mesh, value=data['phase'])
X, Y = mesh.getCellCenters()
ld1=767.671770874
ld2=6586.6684953599997
sd1=7489.1017183100003
sd2=349.289285769
vd1=8.64877867561
vd2=74.207024652200005
vd = vd1 + vd2
ld = ld1 + ld2
phase = (phase > 0.) * phase + (phase > 1.) * ( 1. - phase)
gamma = (density - vd) / (ld - vd)
gamma = (gamma > 0.) * gamma + (gamma > 1.) * ( 1. - gamma)
phis = density1 / density
phil = gamma * (1 - phase)
phiv = (1 - gamma) * (1 - phase)
phis = numerix.reshape(phis, (ny, nx))
phil = numerix.reshape(phil, (ny, nx))
phiv = numerix.reshape(phiv, (ny, nx))
tmps = numerix.zeros((ny, 2 * nx), 'd')
tmpl = numerix.zeros((ny, 2 * nx), 'd')
tmpv = numerix.zeros((ny, 2 * nx), 'd')
tmps[:,:nx] = phis[:,::-1]
tmps[:,nx:] = phis
tmpl[:,:nx] = phil[:,::-1]
tmpl[:,nx:] = phil
tmpv[:,:nx] = phiv[:,::-1]
tmpv[:,nx:] = phiv
phis = tmps
phil = tmpl
phiv = tmpv
## pylab.contour(phis, (0.5,), extent = (-dx * nx * 1e+6, dx * nx * 1e+6, -Ly / 10, dy * ny * 1e+6 - Ly / 10), colors='k')
## pylab.contour(phil, (0.5,), extent = (-dx * nx * 1e+6, dx * nx * 1e+6, -Ly / 10, dy * ny * 1e+6 - Ly / 10), colors='k')
cdict = {'red': ((0.0, 0.0, 0.0),
(0.5, 0.0, 0.0),
(1.0, 0.0, 0.0)),
'green': ((0.0, 0.0, 0.0),
(0.5, 1.0, 1.0),
(1.0, 0.0, 0.0)),
'blue': ((0.0, 0.0, 0.0),
(0.5, 0.0, 0.0),
(1.0, 1.0, 1.0))}
my_cmap = matplotlib.colors.LinearSegmentedColormap('my_colormap',cdict,256)
mm = 1 - phis[::-1]
## pylab.imshow(((1 - phis[::-1]) * 0.5 + phil[::-1] * 0.5) / 2. ,
## extent = (-dx * nx * 1e+6, dx * nx * 1e+6, -Ly / 7, dy * ny * 1e+6 - Ly / 7),
## cmap=my_cmap,
## alpha=1.0)
pylab.imshow(1 - phis[::-1],
extent = (-dx * nx * 1e+6, dx * nx * 1e+6, -Ly / 7, dy * ny * 1e+6 - Ly / 7),
cmap=matplotlib.cm.gray,
alpha=1.0)
pylab.xlim(xmin=-2.5, xmax=2.5)
pylab.xlabel(r'$r$ ($\mu$m)')
pylab.ylabel(r'$h$ ($\mu$m)')
pylab.ylim(ymin=-0.5, ymax=3.5 - 0.5)
pylab.savefig(filename)
time = 0.0
timestep = 0
dataFile = os.path.join('no_tsm', 'data.h5')
data = DictTable(dataFile, 'r')
print data.getLatestIndex()
while timestep <= data.getLatestIndex():
dataFile = os.path.join('no_tsm', 'data.h5.movie')
data = DictTable(dataFile, 'r')
if data[timestep]['elapsedTime'] >= time:
print timestep, time, data[timestep]['elapsedTime']
filename = 'movie/movie' + str(timestep).rjust(7, '0') + '.png'
if not os.path.exists(filename):
makePNG(data[timestep], filename)
time = max(time, data[timestep]['elapsedTime'])
time *= 1.1
timestep += 10
timestep = 800
filename = 'movie/movie' + str(timestep).rjust(7, '0') + '.png'
makePNG(data[timestep], filename)
Raw
parallelCylindricalGrid2D.py
from fipy.meshes.numMesh.cylindricalGrid2D import CylindricalGrid2D
from PyTrilinos import Epetra
from fipy import numerix
class ParallelCylindricalGrid2D(CylindricalGrid2D):
def __init__(self, dx=1., dy=1., nx=None, ny=None, overlap=2, origin=((0.,), (0.,))):
procID = Epetra.PyComm().MyPID()
Nproc = Epetra.PyComm().NumProc()
Ly = dy * ny
self.local_ny = int(ny / Nproc)
local_dy = Ly / Nproc / self.local_ny
local_ly = local_dy * self.local_ny
global_ny = self.local_ny * Nproc
local_origin_y = local_ly * procID + origin[1][0]
## change local values due to overlapping nodes
if Nproc == 1:
overlap = 0
self.overlap = overlap
overlap_local_ny = self.local_ny
if procID > 0:
overlap_local_ny = overlap_local_ny + self.overlap
local_origin_y = local_origin_y - self.overlap * local_dy
if procID < Nproc - 1:
overlap_local_ny = overlap_local_ny + self.overlap
self.globalNumberOfCells = Nproc * self.local_ny * nx
## self.numberOfLocalCellsNoOverlap = self.local_ny * nx
## self.localGlobalIDsNoOverlap = numerix.arange(self.local_ny * nx) + procID * (self.local_ny * nx)
## print 'local_origin_y',local_origin_y
CylindricalGrid2D.__init__(self,
dx=dx,
dy=local_dy,
nx=nx,
ny=overlap_local_ny,
origin=((origin[0][0],), (local_origin_y,)))
def getGlobalCellIDs(self):
procID = Epetra.PyComm().MyPID()
N = self.getNumberOfCells()
localN = self.nx * self.local_ny
shift = self.nx * self.overlap
if procID == 0:
shift = 0
return numerix.arange(N) + localN * procID - shift
def getGlobalLocalCellIDs(self):
procID = Epetra.PyComm().MyPID()
localN = self.nx * self.local_ny
return numerix.arange(localN) + procID * localN
def getGlobalNumberOfCells(self):
return self.globalNumberOfCells
## if procID == 0:
## self.startID = 0
## else:
## self.startID = self.overlap * nx
## if procID == Nproc - 1:
## self.endID = overlap_local_ny * nx - 1
## else:
## self.endID = overlap_local_ny * nx - self.overlap * nx - 1
## def getNumberOfGlobalCells(self):
## return self.numberOfGlobalCells
## def getNumberOfLocalCellsNoOverlap(self):
## return self.numberOfLocalCellsNoOverlap
## def getLocalGlobalIDsNoOverlap(self):
## return self.localGlobalIDsNoOverlap
Raw
radiusVtime.py
from fipy import *
from dictTable import DictTable
from triplePoint import TriplePointContour
import pylab
def getRadiusVtime(dd, tmax, stepsize):
timestep = 0
radiusFile = os.path.join(os.path.join('..', str(dd)), 'radius.gz')
import shutil
datafile = os.path.join(os.path.join(os.path.join('..', str(dd)), 'no_tsm'), 'data.h5')
datafilebkup = datafile + '.bkup'
## shutil.copy(datafile, datafilebkup)
data = DictTable(datafilebkup, 'r')
def getNewVars(timestep):
dx = data[timestep]['dx']
dy = data[timestep]['dy']
nx = data[timestep]['nx']
ny = data[timestep]['ny']
mesh = CylindricalGrid2D(nx=nx, dx=dx, ny=ny, dy=dy)
density1 = CellVariable(mesh=mesh, value=data[timestep]['density1'])
density2 = CellVariable(mesh=mesh, value=data[timestep]['density2'])
phase = CellVariable(mesh=mesh, value=data[timestep]['phase'], name=r'$\phi$')
return density1, density2, phase
density1, density2, phase = getNewVars(0)
delta = 0.1e-6
start = 0.2 * delta
finish = 4. * delta
dis = delta / 2.
Nsteps = data.getLatestIndex() / stepsize
times = []
rads = []
velocities = []
step = 0
tp = TriplePointContour()
if os.path.exists(radiusFile):
times, rads, timestep = dump.read(radiusFile)
else:
times = []
rads = []
timestep = 0
while timestep < data.getLatestIndex() and data[timestep]['elapsedTime'] < tmax:
print 'dir',dd,'timestep',timestep,'t',data[timestep]['elapsedTime']
if len(density1) != len(data[timestep]['density1']):
density1, density2, phase = getNewVars(timestep)
tp = TriplePointContour()
density1.setValue(data[timestep]['density1'])
density2.setValue(data[timestep]['density2'])
phase.setValue(data[timestep]['phase'])
r, y, X, Y = tp.getTriplePointNoContour(phase, density1, density2)
times += [data[timestep]['elapsedTime']]
rads += [r]
timestep += stepsize
dump.write((times, rads, timestep), radiusFile)
return times, rads
pylab.figure()
tmax = 1.e-4
stepsize = 10
Ly = 2.5e-6
R0 = 2. * Ly / 5.
d0 = Ly / 10.0
import numpy
r0 = numpy.sqrt(R0**2 - (Ly / 2. - d0)**2)
## for dd, desc in [(93, '93 central difference N=225'),
## (99, '99 power law N=225'),
## (100, '100 power law N=450'),
## (101, '101 val Leer N=225'),
## (102, '102 power law N=900'),
## (103, '103 van Leer N=450')]:
dd = 124
desc = ''
tt, rr = getRadiusVtime(dd, tmax, stepsize)
print tt, rr
pylab.semilogx(numerix.array(tt) * 1e+6 + 1e-10, (numpy.array(rr) - r0) / R0, label=desc)
##pylab.xlabel(r'$t$ ($\mu$s)')
##pylab.ylabel(r'$\frac{r - r_0}{R_0}$')
##pylab.legend(loc=2)
##pylab.ylabel(r'$\theta_l$ ($\pi$-rad)')
##pylab.legend((r'$\hat{r}=\delta$', r'$\hat{r}=2 \delta$', r'$\hat{r}=3 \delta$'), loc=4)
## ## , r'$\theta_1$', r'$\theta_3$', r'$\theta_4$'), loc=4)
##pylab.text(1.5, 0.34, r'$\theta_l=\theta_l^{\text{equ}}$')
## pylab.text(1.05, 0.2, r'$\delta$')
#ylab.xlim(xmin=1e-4)
pylab.savefig('radiusVtime.png')
Raw
submit_file
Executable = /usr/local/bin/condor_mpi
Arguments = input.py
Universe = parallel
machine_count = 4
Error = $(Cluster).$(Process).$(NODE).err
Output = $(Cluster).$(Process).$(NODE).out
Log = $(Cluster).$(Process).$(NODE).log
Initialdir = /users/wd15/Documents/trunk/reactiveWetting/212
DedicatedScheduler = "DedicatedScheduler@ruth.nist.gov"
environment = "PYTHONPATH=.:/users/wd15/usr/x86_64/4.0/lib/python:/users/wd15/usr/x86_64/4.0/lib/python2.4/site-packages:/users/wd15/usr/x86_64/4.0:/users/wd15/usr/x86_64/4.0/lib/python/gist:/users/wd15/Documents/python/fipy/trunk/utils:/users/wd15/Documents/python/fipy/trunk:/users/wd15/Documents/python/profiler"
getenv = True
Queue
Raw
triplePoint.py
from fipy import *
def getTriplePointField(phase, density, density2, densities1, densities2, phases):
densities1 = numerix.array(densities1)
densities2 = numerix.array(densities2)
phases = numerix.array(phases)
X2s = densities2 / (densities1 + densities2)
densities = densities1 + densities2
density_tj = numerix.sum(densities) / 3.0
phase_tj = numerix.sum(phases) / 3.0
X2_tj = numerix.sum(X2s) / 3.0
phase_denom = max((phases - phase_tj)**2)
X2_denom = max((X2s - X2_tj)**2)
density_denom = max((densities - density_tj)**2)
X2 = density2 / density
return (phase - phase_tj)**2 / phase_denom + \
(X2 - X2_tj)**2 / X2_denom + \
(density - density_tj)**2 / density_denom
## def getArgmin1(density, phase):
## alpha = (density - densities[2]) / (densities[1] - densities[2])
## alpha = (alpha > 0.) * alpha + (alpha > 1.) * ( 1. - alpha)
## print max(phase)
## print min(phase)
## print max(alpha)
## print min(alpha)
## phase_solid = phase
## phase_liquid = alpha * (1. - phase)
## phase_vapor = (1. - alpha) * (1. - phase)
## print max(phase_solid + phase_liquid + phase_vapor)
## print min(phase_solid + phase_liquid + phase_vapor)
## raw_input()
## third = 1. / 3.
## field = (phase_solid - third)**2 + (phase_liquid - third)**2 + (phase_vapor - third)**2
## argmin = numerix.argmin(field)
## return argmin
def getArgmin(density1, density2, phase, densities1, densities2, phases):
density1_tj = numerix.sum(densities1) / 3.0
density2_tj = numerix.sum(densities2) / 3.0
phase_tj = numerix.sum(phases) / 3.0
phase_denom = max((phases - phase_tj)**2)
density1_denom = max((densities1 - density1_tj)**2)
density2_denom = max((densities2 - density2_tj)**2)
field = (phase - phase_tj)**2 / phase_denom + \
(density1 - density1_tj)**2 / density1_denom + \
(density2 - density2_tj)**2 / density2_denom
return numerix.argmin(field)
def getTriplePoint(phase, density1, density2, densities1, densities2, phases):
densities1 = numerix.array(densities1)
densities2 = numerix.array(densities2)
phases = numerix.array(phases)
x, y = density2.getMesh().getCellCenters()[:,getArgmin(density1, density2, phase, densities1, densities2, phases)]
## print x, y
## raw_input('course')
dx = density2.getMesh().dx
dy = density2.getMesh().dy
refine = 10
subMesh = Grid2D(nx=refine, ny=refine, dx=2 * dx / refine, dy=2 * dy / refine)
origin = numerix.zeros((2, len(subMesh.getFaceCenters()[0])), 'd')
origin[0] = x - dx
origin[1] = y - dy
centers = subMesh.getFaceCenters() + origin
density1 = density1(centers, order=1)
density2 = density2(centers, order=1)
phase = phase(centers, order=1)
return centers[:,getArgmin(density1, density2, phase, densities1, densities2, phases)]
##def __getData(var):
## from fipy.tools.numerix import array, reshape
## return reshape(array(var), var.getMesh().getShape()[::-1])
class TriplePointContour:
def __getData(self, var, mesh):
from fipy.tools.numerix import array, reshape
return reshape(var, mesh.getShape()[::-1])
def getTriplePointContour(self, phase, density1, density2):
import pylab
Lx, Ly = phase.getMesh().getPhysicalShape()
nx, ny = phase.getMesh().getShape()
dx = Lx / nx
dy = Ly / ny
shift = dx / 10.0
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + ((shift,),(0,))
density = density1 + density2
density = CellVariable(mesh=mesh, value=density)
phase = CellVariable(mesh=mesh, value=phase)
density1 = CellVariable(mesh=mesh, value=density1)
density2 = CellVariable(mesh=mesh, value=density2)
N = 2000
newdx = Lx / N
newMesh = Grid2D(nx=N, ny=N, dx=newdx, dy=newdx) + ((shift,),(0,))
newPhase = numerix.zeros(newMesh.getNumberOfCells(), 'd')
newDensity = numerix.zeros(newMesh.getNumberOfCells(), 'd')
if not hasattr(self, 'nearestCellIDs'):
self.cc = newMesh.getCellCenters()
self.nearestCellIDs = mesh._getNearestCellID(self.cc)
newPhase[:] = numerix.array(phase(self.cc, order=1, nearestCellIDs=self.nearestCellIDs))
newDensity[:] = numerix.array(density(self.cc, order=1, nearestCellIDs=self.nearestCellIDs))
ld1=767.671770874
ld2=6586.6684953599997
sd1=7489.1017183100003
sd2=349.289285769
vd1=8.64877867561
vd2=74.207024652200005
vd = vd1 + vd2
ld = ld1 + ld2
gamma = (newDensity - vd) / (ld - vd)
gamma = (gamma > 0.) * gamma + (gamma > 1.) * ( 1. - gamma)
import pylab
pylab.figure()
newPhase = (newPhase > 0.) * newPhase + (newPhase > 1.) * ( 1. - newPhase)
phi1 = newPhase
phi2 = gamma * (1 - newPhase)
phi3 = (1 - gamma) * (1 - newPhase)
extent = (shift, Lx + shift, 0, Ly)
third = 1. / 3.
field = (phi1 - third)**2 + (phi2 - third)**2 + (phi3 - third)**2
minValue = field[numerix.argmin(field)]
c4 = pylab.contour(self.__getData(field, newMesh),
(minValue + 1e-2,),
extent=extent,
alpha=1.0,
colors='black')
c5 = pylab.contour(self.__getData(newPhase, newMesh),
(third,),
extent=extent,
alpha=1.0,
colors='black')
vv = numerix.array(c4.collections[0].get_verts(), 'd')
vv = vv[numerix.argmin(vv[:,0])]
bottom = numerix.array(c5.collections[0].get_verts(), 'd')
bottom = bottom[numerix.argmin(bottom[:,0])]
return vv[0], vv[1], bottom[0], bottom[1]
def getTriplePointNoContour(self, phase, density1, density2):
import pylab
Lx, Ly = phase.getMesh().getPhysicalShape()
nx, ny = phase.getMesh().getShape()
dx = Lx / nx
dy = Ly / ny
shift = dx / 10.0
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + ((shift,),(0,))
density = density1 + density2
density = CellVariable(mesh=mesh, value=density)
phase = CellVariable(mesh=mesh, value=phase)
density1 = CellVariable(mesh=mesh, value=density1)
density2 = CellVariable(mesh=mesh, value=density2)
N = 2000
newdx = Lx / N
newMesh = Grid2D(nx=N, ny=N, dx=newdx, dy=newdx) + ((shift,),(0,))
newPhase = numerix.zeros(newMesh.getNumberOfCells(), 'd')
newDensity = numerix.zeros(newMesh.getNumberOfCells(), 'd')
if not hasattr(self, 'nearestCellIDs'):
self.cc = newMesh.getCellCenters()
self.nearestCellIDs = mesh._getNearestCellID(self.cc)
newPhase[:] = numerix.array(phase(self.cc, order=1, nearestCellIDs=self.nearestCellIDs))
newDensity[:] = numerix.array(density(self.cc, order=1, nearestCellIDs=self.nearestCellIDs))
ld1=767.671770874
ld2=6586.6684953599997
sd1=7489.1017183100003
sd2=349.289285769
vd1=8.64877867561
vd2=74.207024652200005
vd = vd1 + vd2
ld = ld1 + ld2
gamma = (newDensity - vd) / (ld - vd)
gamma = (gamma > 0.) * gamma + (gamma > 1.) * ( 1. - gamma)
import pylab
pylab.figure()
newPhase = (newPhase > 0.) * newPhase + (newPhase > 1.) * ( 1. - newPhase)
phi1 = newPhase
phi2 = gamma * (1 - newPhase)
phi3 = (1 - gamma) * (1 - newPhase)
extent = (shift, Lx + shift, 0, Ly)
half = 3. / 5.
quarter = 1. / 5.
field = (phi1 - quarter)**2 + (phi2 - half)**2 + (phi3 - quarter)**2
minValue = field[numerix.argmin(field)]
if not hasattr(self, 'newMeshCellCenters'):
self.newMeshCellCenters = newMesh.getCellCenters()
x,y = self.newMeshCellCenters[:,numerix.argmin(field)]
self.liquidContour = pylab.contour(self.__getData(phi2, newMesh),
(half,), extent=extent, oalpha=1.0, colors='black')
self.solidContour = pylab.contour(self.__getData(phi1, newMesh),
(quarter,), extent=extent, alpha=1.0, colors='black')
self.vaporContour = pylab.contour(self.__getData(phi3, newMesh),
(quarter,), extent=extent, alpha=1.0, colors='black')
## bottom = numerix.array(self.solidContour.collections[0].get_verts(), 'd')
## bottom = bottom[numerix.argmin(bottom[:,0])]
self.triplePoint = (x, y)
self.phi2 = phi2
self.newMesh = newMesh
## return x, y, bottom[0], bottom[1]
return x, y, None, None
def getPoints(self, radius):
vvs = []
xlow = []
ylow = []
xhigh = []
yhigh = []
for cc in self.vaporContour, self.solidContour, self.liquidContour:
vv = numerix.array(cc.collections[0].get_verts(), 'd')
Rfunc = abs(radius - numerix.sqrt((self.triplePoint[0] - vv[:,0])**2 + (self.triplePoint[1] - vv[:,1])**2))
as = numerix.argsort(Rfunc)
p0 = [vv[as[0],0], vv[as[0],1]]
p1 = [vv[as[1],0], vv[as[1],1]]
nx = numerix
if nx.sqrt(nx.sum((nx.array(p0)- nx.array(p1))**2)) < (radius / 10.):
p1 = [vv[as[2],0], vv[as[2],1]]
## print 'got here1'
## if nx.sqrt(nx.sum((nx.array(p0)- nx.array(p1))**2)) < (radius / 10.):
## p1 = [vv[as[3],0], vv[as[3],1]]
## print 'got here2'
vvs += [[p0, p1]]
if p0[0] < p1[0]:
xlow += [0]
xhigh += [1]
else:
xlow += [1]
xhigh += [0]
if p0[1] < p1[1]:
ylow += [0]
yhigh += [1]
else:
ylow += [1]
yhigh += [0]
import numpy
solidVaporPoint = (numpy.array(vvs[0][ylow[0]]) + numpy.array(vvs[1][xhigh[1]])) / 2
liquidVaporPoint = (numpy.array(vvs[0][yhigh[0]]) + numpy.array(vvs[2][yhigh[2]])) / 2
solidLiquidPoint = (numpy.array(vvs[1][xlow[1]]) + numpy.array(vvs[2][xlow[2]])) / 2
## print 'solidVaporPoint',solidVaporPoint
## raw_input("stop")
return [solidVaporPoint, liquidVaporPoint, solidLiquidPoint]
## return vvs
def getAngles(self, radius, phase, density1, density2, newTriplePoint=False):
if not hasattr(self,'_tp'):
xtp, ytp, bx, by = self.getTriplePointNoContour(phase, density1, density2)
self._tp = numerix.array((xtp, ytp), 'd')
sv, lv, sl = self.getPoints(radius)
def calcAngle(p1, p2):
nx = numerix
c = nx.sqrt(nx.sum((p1 - p2)**2))
a = nx.sqrt(nx.sum((p1 - self._tp)**2))
b = nx.sqrt(nx.sum((p2 - self._tp)**2))
return nx.arccos((a**2 + b**2 - c**2) / (2 * a * b))
## print 'radius',radius
## print 'xtp, ytp',xtp, ytp
## print 'sv',sv
## print 'lv',lv
## print 'sl',sl
## print 'angles',numerix.array((calcAngle(sv, lv),
## calcAngle(lv, sl),
## calcAngle(sl, sv))) * 180. / numerix.pi
## raw_input("stopped")
return (calcAngle(sv, lv),
calcAngle(lv, sl),
calcAngle(sl, sv))
def getAngles1(self, radius, dis, phase, density1, density2, newTriplePoint=True):
## if not hasattr(self,'_tp'):
xtp, ytp, bx, by = self.getTriplePointNoContour(phase, density1, density2)
self._tp = numerix.array((xtp, ytp), 'd')
sv1, lv1, sl1 = self.getPoints(radius)
sv2, lv2, sl2 = self.getPoints(radius + dis)
def calcAngle(p1, p2, tp):
nx = numerix
c = nx.sqrt(nx.sum((p1 - p2)**2))
a = nx.sqrt(nx.sum((p1 - tp)**2))
b = nx.sqrt(nx.sum((p2 - tp)**2))
return nx.arccos((a**2 + b**2 - c**2) / (2 * a * b))
np = lv1.copy()
np[0] = np[0] - dis
theta1 = calcAngle(np, lv2, lv1)
theta2 = numerix.pi - theta1
## print np, lv2, lv1
## print 'radius',radius
## print 'dis',dis
## print 'theta1',theta1
## raw_input("stop")
np = sl1.copy()
np[0] = np[0] - dis
theta3 = calcAngle(np, sl2, sl1)
if sl2[1] > sl1[1]:
theta3 = - theta3
np = sv1.copy()
np[0] = np[0] + dis
theta4 = calcAngle(np, sv2, sv1)
if sv2[1] > sv1[1]:
theta4 = -theta4
theta5 = numerix.pi - theta3 - theta4
## print 'radius',radius
## print 'xtp, ytp',xtp, ytp
## print 'sv',sv
## print 'lv',lv
## print 'sl',sl
## print 'angles',numerix.array((calcAngle(sv, lv),
## calcAngle(lv, sl),
## calcAngle(sl, sv))) * 180. / numerix.pi
## raw_input("stopped")
return theta2 + theta4, theta1 + theta3, theta5, theta1, theta3, theta4
Raw
triplePointVelocityVtime.py
from dictTable import DictTable
from matplotlib import rc, rcParams
import matplotlib
matplotlib.use('Agg')
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rc('text', usetex=True)
import pylab
pylab.ioff()
pylab.hold(True)
from fipy import *
import numpy
import scipy.interpolate
dd = 119
radiusFile = 'radius.gz'
times, rads, timestep = dump.read(radiusFile)
times = numpy.array(times)
N = 10000
times2 = 10**(-10. + numpy.arange(N) * (numpy.log10(times[-1]) + 11.) / (N - 1))
times2 = times2[times2 < times[-1]]
weights = numpy.ones(len(times))
##weights[1:] = numpy.log(times[1:]) / numpy.log(times[1])
print weights[0]
print weights[-1]
##weights[weights < 1e-20] = 1e-20
##weights[weights > 1e+20] = 1e+20
rads2 = scipy.interpolate.splev(times2,
scipy.interpolate.splrep(times,
rads,
k=3,
w=weights,
s=4e-15))
pylab.semilogx(times2, rads2, times, rads)
pylab.savefig('tmp.png')
v2 = scipy.interpolate.splev(times2,
scipy.interpolate.splrep(times,
rads,
k=3,
w=weights,
s=4e-15),
der=1)
pylab.figure()
Ly = 2.5e-6
R0 = 2. * Ly / 5.
viscosity = 2e-3
density = 6586.6684953599997
gamma = (19.114 + 18.904) / 2.
oh = viscosity / numerix.sqrt(density * gamma * R0)
print 'oh',oh
sliplength = 0.05e-6
lamb = sliplength / R0
print 'lamb',lamb
T = numerix.sqrt(density * R0**3 / gamma) * (density * gamma * R0 / viscosity**2)**(1. / 8.)
Udim = R0 / T
t4 = times2 / T
v4 = density * v2 * sliplength / viscosity
re = v2 * density * sliplength / viscosity
RE = sliplength / R0
pylab.loglog(t4, re, 'b', label=r'this work')
t5 = t4[(t4 > 0.02) & (t4 < 1)]
pylab.loglog(t5, 10.0 * t5**-0.5 / 25.0, 'k--', label=r'$t^{-0.5}$')
##t7 = t4[(t4 > 1) & (t4 < 10)]
##pylab.loglog(t7, 10.0 * t7**-0.8 / 25.0, 'k-.', label=r'$t^{-0.8}$')
##t6 = t4[t4 > 10]
##pylab.loglog(t6, 10.0 * t6**-2 / 1.5, 'k:', label=r'$t^{-2}$')
exp = numpy.loadtxt('download.csv', skiprows=1, delimiter=',')
exptime = 10**exp[:,0]
expU = 10**exp[:,1] / 0.01 * 0.00707
pylab.loglog(exptime, expU, 'r', label=r'Ding \it{et. al.}')
pylab.xlabel(r'$t / T$')
pylab.ylabel(r'$\lambda Re$')
pylab.loglog((5., 50.), (RE, RE), 'k:', label=r'$Re=1$')
pylab.legend()
pylab.savefig('triplePointVelocityVtime.pdf')
Raw
unsegregatedEquation.py
#!/usr/bin/env python
from fipy import *
from PyTrilinos import Epetra
class UnsegregatedEquation:
def __init__(self, equations):
self.equations = equations
def buildMatrix(self, vars, dt, boundaryConditions):
N = vars[0].getMesh().getNumberOfCells()
M = len(vars)
solutionVector = numerix.zeros(M * N, 'd')
matrixdiag = numerix.zeros(M * N, 'd')
RHSvector = numerix.zeros(M * N, 'd')
from fipy.tools.pysparseMatrix import _PysparseMatrix
matrix = _PysparseMatrix(size=M * N)
for i in range(M):
solutionVector[i * N:(i + 1) * N] = vars[i]
for j in ((numerix.arange(M) - (M - i)) % M):
if self.equations[i][j] is not None:
s, tmpMatrix, tmpRHSvector = self.equations[i][j]._prepareLinearSystem(vars[j], None, boundaryConditions[i][j], dt)
if i == j:
diag = tmpMatrix.takeDiagonal()
jacobi = tmpMatrix.__class__(size=N)
JIDS = numerix.arange(N)
jacobi.addAt(1 / diag, JIDS, JIDS)
## jacobi.putDiagonal(1 / diag) ## memory leak
matrixdiag[i * N:(i + 1) * N] = diag
tmpRHSvector = jacobi * numerix.array(tmpRHSvector)
tmpMatrix = jacobi * tmpMatrix
RHSvector[i * N:(i + 1) * N] += tmpRHSvector
del tmpRHSvector
keys = tmpMatrix.matrix.keys()
if len(keys) > 0:
ids = numerix.zeros((len(keys[0]), 2), 'l')
ids[:,0] = keys[0]
ids[:,1] = keys[1]
ids[:,0] = ids[:,0] + i * N
ids[:,1] = ids[:,1] + j * N
matrix.addAt(tmpMatrix.matrix.values(), ids[:,0], ids[:,1])
del ids
del keys
return matrix, solutionVector, RHSvector, matrixdiag
def calcIDs(self, IDs, M, mesh):
N = len(IDs)
IDs = numerix.resize(IDs, (M, N))
tmp = numerix.zeros((M, N), 'l')
for m in range(M):
tmp[m] = m * mesh.getGlobalNumberOfCells()
return (IDs + tmp).flatten()
def calcLocal(self, vec, M, N, mesh, comm):
Nproc = comm.NumProc()
if Nproc == 1:
return vec
else:
procID = comm.MyPID()
overlap = mesh.overlap * mesh.nx
tmp = numerix.reshape(vec, (M, N))
if procID == 0:
return tmp[:, :-overlap].flatten()
elif procID == Nproc - 1:
return tmp[:, overlap:].flatten()
else:
return tmp[:, overlap:-overlap].flatten()
def calcMatrix(self, matrix, M, N, mesh, map, comm):
trilinosMatrix = Epetra.CrsMatrix(Epetra.Copy, map, -1)
Nproc = comm.NumProc()
procID = comm.MyPID()
totalCells = Nproc * mesh.local_ny * mesh.nx
for i in range(M):
for j in range(M):
NOC = mesh.overlap * mesh.nx
i0 = i * N
i1 =(i + 1) * N
j0 = j * N
j1 = (j + 1) * N
if procID == 0:
A = matrix.matrix[i0:i1 - NOC,j0:j1]
elif procID == Nproc - 1:
A = matrix.matrix[i0 + NOC:i1,j0:j1]
else:
A = matrix.matrix[i0 + NOC:i1 - NOC,j0:j1]
keys = A.keys()
## print type(keys[0][0])
## raw_input("stop")
IDs = numerix.zeros((len(keys[0]), 2), 'l')
IDs[:,0] = numerix.array(keys[0], 'l')
IDs[:,1] = numerix.array(keys[1], 'l')
IDs[:,0] += (i * totalCells + procID * mesh.local_ny * mesh.nx)
if procID == 0:
IDs[:,1] += (j * totalCells + procID * mesh.local_ny * mesh.nx)
else:
IDs[:,1] += (j * totalCells + procID * mesh.local_ny * mesh.nx - mesh.nx * mesh.overlap)
trilinosMatrix.InsertGlobalValues(IDs[:,0].astype('int32'), IDs[:,1].astype('int32'), A.values())
return trilinosMatrix
def sweep(self, vars=vars, dt=1.0, boundaryConditions=(), relaxation=0.5):
M = len(vars)
matrix, solutionVector, RHSvector, matrixdiag = self.buildMatrix(vars, dt, boundaryConditions)
mesh = vars[0].getMesh()
N = mesh.getNumberOfCells()
IDs = self.calcIDs(mesh.getGlobalCellIDs(), M, mesh)
comm = Epetra.PyComm()
global_map = Epetra.Map(-1, list(IDs), 0, comm)
globalVector = Epetra.Vector(global_map)
globalVector[:] = solutionVector
local_IDs = self.calcIDs(mesh.getGlobalLocalCellIDs(), M, mesh)
local_map = Epetra.Map(-1, list(local_IDs), 0, comm)
local_solutionVector = self.calcLocal(solutionVector, M, N, mesh, comm)
local_solutionVector = Epetra.Vector(local_map, local_solutionVector)
local_RHSvector = self.calcLocal(RHSvector, M, N, mesh, comm)
local_RHSvector = Epetra.Vector(local_map, local_RHSvector)
## build matrix
trilinosMatrix = self.calcMatrix(matrix, M, N, mesh, local_map, comm)
trilinosMatrix.FillComplete()
trilinosMatrix.OptimizeStorage()
## calculate residual
residual = Epetra.Vector(local_map)
tmp = Epetra.Vector(local_map)
trilinosMatrix.Multiply(False, local_solutionVector, tmp)
residual[:] = tmp - local_RHSvector
residual = residual.Norm2()
procID = comm.MyPID()
## solve matrix
from PyTrilinos import AztecOO
solver = AztecOO.AztecOO(trilinosMatrix, local_solutionVector, local_RHSvector)
solver.SetAztecOption(AztecOO.AZ_solver, AztecOO.AZ_gmres)
solver.SetAztecOption(AztecOO.AZ_output, AztecOO.AZ_none)
solver.SetAztecOption(AztecOO.AZ_precond, AztecOO.AZ_Jacobi)
## solver.SetAztecOption(AztecOO.AZ_precond, AztecOO.AZ_none)
solver.SetAztecOption(AztecOO.AZ_poly_ord, 2)
solver.SetAztecOption(AztecOO.AZ_kspace, 30)
convergence = solver.Iterate(2000, 1e-4)
## redistribute
globalVector.Import(local_solutionVector, Epetra.Import(global_map, local_map), Epetra.Insert)
for i in range(M):
vars[i].setValue(relaxation * globalVector[i * N:(i + 1) * N] + (1 - relaxation) * vars[i])
return residual, matrixdiag, convergence
Raw
velocity.py
from fipy import *
from dictTable import DictTable
from triplePoint import TriplePointContour
from matplotlib import rc, rcParams
import matplotlib
matplotlib.use('Agg')
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rc('text', usetex=True)
import pylab
pylab.ioff()
pylab.hold(True)
pylab.figure()
timestep = 1100
dataFile = os.path.join('no_tsm', 'data.h5')
data = DictTable(dataFile, 'r')
dx = data[timestep]['dx']
dy = data[timestep]['dy']
nx = data[timestep]['nx']
ny = data[timestep]['ny']
mesh = CylindricalGrid2D(nx=nx, dx=dx, ny=ny, dy=dy)
Xvelocity = CellVariable(mesh=mesh, value=data[timestep]['Xvelocity'])
Yvelocity = CellVariable(mesh=mesh, value=data[timestep]['Yvelocity'])
density1 = CellVariable(mesh=mesh, value=data[timestep]['density1'])
density2 = CellVariable(mesh=mesh, value=data[timestep]['density2'])
phase = CellVariable(mesh=mesh, value=data[timestep]['phase'])
density = density1 + density2
Xvelocity = Xvelocity * density
Yvelocity = Yvelocity * density
X, Y = mesh.getCellCenters()
##timestep=1000
Xrange = (0.6e-6, 1.0e-6)
Yrange = (0.1e-6, 0.5e-6)
Xmask = (X > Xrange[0]) & (X < Xrange[1])
Ymask = (Y > Yrange[0]) & (Y < Yrange[1])
mask = Xmask & Ymask
pylab.quiver(X[mask], Y[mask], Xvelocity[mask], Yvelocity[mask], scale=0.5e+6)
## N = 2000
## newdx = Lx / N
## shift = dx / 10.0
## newMesh = Grid2D(nx=N, ny=N, dx=newdx, dy=newdx) + ((shift,),(0,))
## newPhase = numerix.zeros(newMesh.getNumberOfCells(), 'd')
## newDensity = numerix.zeros(newMesh.getNumberOfCells(), 'd')
## cc = newMesh.getCellCenters()
## nearestCellIDs = mesh._getNearestCellID(self.cc)
## phase = numerix.array(phase(self.cc, order=1, nearestCellIDs=self.nearestCellIDs))
## density = numerix.array(density(self.cc, order=1, nearestCellIDs=self.nearestCellIDs))@
ld1=767.671770874
ld2=6586.6684953599997
sd1=7489.1017183100003
sd2=349.289285769
vd1=8.64877867561
vd2=74.207024652200005
vd = vd1 + vd2
ld = ld1 + ld2
phase = (phase > 0.) * phase + (phase > 1.) * ( 1. - phase)
gamma = (density - vd) / (ld - vd)
gamma = (gamma > 0.) * gamma + (gamma > 1.) * ( 1. - gamma)
phis = phase
phil = gamma * (1 - phase)
phiv = (1 - gamma) * (1 - phase)
##nx = N
##ny = N
##nx = int((Xrange[1] - Xrange[0]) / dx + 1e-10)
##ny = int((Yrange[1] - Yrange[0]) / dy + 1e-10)
##phis = numerix.reshape(phis[mask], (nx, ny))
##phil = numerix.reshape(phil[mask], (nx, ny))
##pylab.contour(phis, (0.5,), extent = (Xrange[0], Xrange[1], Yrange[0], Yrange[1]))
##pylab.contour(phil, (0.1,), extent = (Xrange[0], Xrange[1], Yrange[0], Yrange[1]))
##timestep=10000
##mask = (X > 1.0e-6) & (X < 2.0e-6) & (Y > 0.0e-6) & (Y < 1.0e-6)
##pylab.quiver(X[mask], Y[mask], Xvelocity[mask], Yvelocity[mask], scale=0.5e+1)
pylab.savefig('velocity.png')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment