ver = 'NSThermN_v1.10'    # FiPy script version
nparts = 100              # number of nanoparticles in simulation
nsteps = 500              # number of time steps
aintv = 5                 # time interval between data averages
pintv = 5                 # time interval between data plot output
#                         #  i.e., (tstep+1) modulo pintv
 
## FiPy script
## to elucidate heat transfer from laser-heated gold-coated 
## silica-core nanoparticles in a homogeneous medium at 37 C
#
## all lengths in units of microns
## all temperatures in Celsius
#
## The gold coating may have a pie-shaped wedge hole 
## along the x'-axis: wedge angle = 0.8 radians; 
## arc length at 65 nm = 52 nm. The x'-axis is at 
## a random angle with respect to the horizontal

from fipy import *

# simulation size
# express all lengths in units of microns
nm = 1e-3
micron = 1
Lx = 5000*nm
Ly = 5000*nm

# simulation resolution
dx = 5*nm		# nx = 1000
dy = 5*nm		# ny = 1000

#dx = 10*nm		# nx = 500
#dy = 10*nm		# ny = 500

#dx = 20*nm		# nx = 250
#dy = 20*nm		# ny = 250

nx = Lx / dx
ny = Ly / dy

mesh = Grid2D(dx = dx, dy=dy, nx = nx, ny=ny)

# create the temperature field
temp = CellVariable(mesh = mesh, name = 'temperature', value = 0.0, hasOld = 1.0)

# create phase fields
phi = CellVariable(mesh = mesh, name = 'phase', value = 0.0, hasOld = 0.0)
# note: the decimal point is critical in the "value" value, or phi is 
# interpreted as an integer with incorrect results for what is wanted."

# physical properties for water medium
rhom = 1.0e-6     # density [micrograms /(micron)^3]
Cm = 4181.3       # heat capacity [nJ/microgram/K]
rhoCm = rhom * Cm
km = 0.581e-3     # thermal conductivity [mW/micron/K]

# physical properties for gold (Au) shell
rhos = 19.3e-6    # density [micrograms /(micron)^3]
Cs = 129.1        # heat capacity [nJ/microgram/K]
rhoCs = rhos * Cs
ks = 320.0e-3     # thermal conductivity [mW/micron/K]

# physical properties for silica particle cores
rhoc = 2.2e-6     # density [micrograms /(micron)^3]
Cc = 703.0        # heat capacity [nJ/microgram/K]
rhoCc = rhoc * Cc
kc = 1.4e-3       # thermal conductivity [mW/micron/K]

# create a source field
source = CellVariable(mesh = mesh, name = 'heat_source', value = 0.0, hasOld = 1.0)

# heat source rate
express all power in units of milli Watts
source_rate = 1.0    # 1.0 mW / micron^3 = 1.0 nJ / microsecond / micron^3


x, y = mesh.getCellCenters()

core_yes=0
shell_yes=0
ts = 15.0*nm     # shell thickness
Rc = 50.0*nm     # particle core radius
Rp = Rc + ts     # particle radius
zetad = 0.4      # wedge half-angle of coating defect

# define particle region: part_yes
part_yes=0
pxmin = 0.2*Lx
pxmax = 0.8*Lx
pymin = 0.2*Ly
pymax = 0.8*Ly

part_yes = ((x>pxmin)&(x<pxmax)&(y>pymin)&(y<pymax))|(part_yes)

cxmin = pxmin + Rp
cxrgn = pxmax - Rp - cxmin
cymin = pymin + Rp
cyrgn = pymax - Rp - cymin

# place particles
random.seed(94763479)    # set random number seed #0
#random.seed(76953371)    # set random number seed #1
#random.seed(20536404)    # set random number seed #2
#random.seed(72970159)    # set random number seed #3
#random.seed(90481135)    # set random number seed #4

#random.seed(42789210)    # set random number seed #5
#random.seed(58330515)    # set random number seed #6
#random.seed(68125755)    # set random number seed #7
#random.seed(72562932)    # set random number seed #8
#random.seed(34404799)    # set random number seed #9

# place nanoparticles in particle region
for i in range(nparts):
    vx = cxmin + cxrgn*random.random()
    vy = cymin + cyrgn*random.random()

##   algorithm to place circular particles
#    tx = x - vx
#    ty = y - vy
#    core_yes = (tx*tx+ty*ty<=Rp*Rp)|(core_yes)

#   algorithm to place circular particles with core and shell
    theta=numerix.pi*2*random.random()
    c=numerix.cos(theta)
    s=numerix.sin(theta)
#   coordinates for a randomly rotated particle
    tx=c*(x-vx)+s*(y-vy)
    ty= -s*(x-vx)+c*(y-vy)
    core_yes=(tx*tx+ty*ty<=Rc*Rc)|(core_yes)

#    shell without a shell defect
#    shell_yes = (tx*tx+ty*ty>Rc*Rc)&(tx*tx+ty*ty<=Rp*Rp)|(shell_yes)

#   shell with a single shell defect
#   cylindrial angle for point (tx, ty)
    zeta = numerix.arctan2(ty, tx)
    shell_yes=(tx*tx+ty*ty>Rc*Rc)&(tx*tx+ty*ty<=Rp*Rp)&((zeta>zetad)|(zeta<-zetad))|(shell_yes)

# initial conditions for particles
temp.setValue(37.0)
temp.setValue(37.0,where=shell_yes)
temp.setValue(37.0,where=core_yes)

phi.setValue(0.50)
phi.setValue(0.75,where=part_yes)
phi.setValue(1.25,where=shell_yes&~core_yes)
phi.setValue(1.75,where=core_yes)

# compute area fractions
# area of pixel: dx * dy
sim_area = (phi >= 0).getCellVolumeAverage()    # simulation area in pixels
part_rgn = (phi > 0.7).getCellVolumeAverage()   # area of particle region in pixels
parts = (phi > 1.2).getCellVolumeAverage()      # area of particles in pixels
cores = (phi > 1.7).getCellVolumeAverage()      # area of cores in pixels
shells = parts - cores                          # area of shells in pixels 
part_f = parts / part_rgn
core_f = cores / part_rgn
shell_f = shells / part_rgn
rgn_f = 1.0 - part_f

source.setValue(0.0)
source.setValue(source_rate,where=shell_yes&~core_yes)
source.setValue(0.0,where=core_yes)

# show mesh
#viewer_m = make(vars = (phi,),limits={'datamin': 0.0, 'datamax': 2.5}, palette='jet.gp')
viewer_m = make(vars = (phi,),limits={'datamin': 0.0, 'datamax': 2.0})
viewer_m.plot(filename='mesh.png')

# dump mesh
dump.write(phi,filename='phase.gz')   # dump (write) mesh to a file for later reading
#dump.read(phi,filename='phase.gz')   # read phase (phi) and analyze and/or plot

# set properties
conduc = CellVariable(mesh=mesh, name = 'conductivity')
conduc.setValue(km)
conduc.setValue(ks, where=shell_yes)
conduc.setValue(kc, where=core_yes)

rhoC = CellVariable(mesh=mesh, name = 'density_heatcapacity')
rhoC.setValue(rhoCm)
rhoC.setValue(rhoCs, where=shell_yes)
rhoC.setValue(rhoCc, where=core_yes)

tempeq = TransientTerm(coeff=rhoC) - ImplicitDiffusionTerm(coeff=conduc) - source

##solver = LinearLUSolver()
solver = LinearPCGSolver()
##solver = LinearCGSSolver()
#
# ##########################
# #  Begin full exclusion  #
# ##########################

import datetime

#raw_input()

maxsweep=100
# 
# express all times in units of microseconds
#dt = 2.0e-6   #  2.0 ps
#dt = 0.002    #  2.0 ns
dt = 0.010    # 10.0 ns
#dt = 0.020    # 20.0 ns

# output log
WRT = open(ver+'.log', "a")
print >> WRT, 'version: ', ver
print >> WRT, 'file written', datetime.datetime.now()
print >> WRT, 'particle frac = ', part_f, '   particle region frac. = ', rgn_f
print >> WRT, '    core frac = ', core_f, '   shell frac. = ', shell_f
print >> WRT, 'no. particles = ', nparts, '\n'
WRT.close()

# viewer for temperature
viewer_s = make(vars = (temp,),limits={'datamin': 35.0, 'datamax': 65.0})
viewer_s.plot()
#from fipy.viewers.tsvViewer import TSVViewer

for tstep in range(nsteps):
    temp.updateOld()
#    source.updateOld()
    tempeq.solve(var=temp,dt=dt,solver=solver)

    time = (tstep+1)*dt

    if (tstep+1)%aintv==0:
        avetemp = temp.getCellVolumeAverage() / sim_area
        partemp = ((phi>0.7)*temp).getCellVolumeAverage() / part_rgn

        WRT = open(ver+'.log', "a")
        print >> WRT, 'step',tstep,'time=',time,'max_temp=',max(temp), \
                  'ave(part_temp)=',partemp,'ave(temp)=',avetemp,'min_temp=',min(temp)
        WRT.close()

    if (tstep+1)%pintv==0:
        viewer_m.plot()
        viewer_s.plot(filename='step%03d.png' % tstep)

#        TSVViewer(vars=(temp, temp.getGrad())).plot(filename='step%03d.dat.gz' % tstep)

#        dump.write(temp,filename='step%03d.gz' % tstep)   # to write temp to a file for later reading
#        dump.read(temp,filename='step%03d.gz' % tstep)   # to read temp later and analyze and/or plot

WRT = open(ver+'.log', "a")
print >> WRT, '\n', 'file written', datetime.datetime.now()
WRT.close()