#!/usr/bin/env python## This script was derived from## 'examples/phase/impingement/mesh20x20/input.py'# In the following examples, we solve the same set of equations as in::##     $ examples/phase/impingement/mesh40x1/input.py## with different initial conditions and a 2D mesh:## .. raw:: latex##    \IndexModule{parser}## ..#from fipy.tools.parser import parse
#numberOfElements = parse('--numberOfElements', action = 'store',
                         type = 'int', default = 40000)
numberOfSteps = parse('--numberOfSteps', action = 'store',
                      type = 'int', default = 100)
## .. raw:: latex##    \IndexModule{numerix}#    \IndexFunction{sqrt}#    \IndexClass{Grid2D}## ..#steps = numberOfSteps
from fipy.tools.numerix import sqrt
N = int(sqrt(numberOfElements))
L = 2.5 * N / 100.
dL = L / N
dx = 0.02
dy = 0.02
from fipy.meshes.grid2D import Grid2D
mesh = Grid2D(dx=dx, dy=dx, nx=N, ny=N)
## The initial conditions are given by## .. raw:: latex##     $ \phi = 1 $ and##     $$ \theta = \begin{cases}#     \frac{2 \pi}{3} & \text{for $x^2 - y^2 < L / 2$,} \\#     \frac{-2 \pi}{3} & \text{for $(x-L)^2 - y^2 < L / 2$,} \\#     \frac{-2 \pi}{3}+0.3 & \text{for $x^2 - (y-L)^2 < L / 2$,} \\#     \frac{2 \pi}{3} & \text{for $(x-L)^2 - (y-L)^2 < L / 2$.}#     \end{cases} $$## This defines four solid regions with different# orientations. Solidification occurs and then boundary wetting occurs# where the orientation varies.## The parameters for this example are#timeStepDuration = 0.02
phaseTransientCoeff = 0.1
thetaSmallValue = 1e-6
beta = 1e5
mu = 1e3
thetaTransientCoeff = 0.01
gamma= 1e3
epsilon = 0.008
s = 0.01
alpha = 0.015
initial_temperature= 10.0
## The system is held isothermal at#
## and is initialized to liquid everywhere## .. raw:: latex##    \IndexClass{CellVariable}## ..#from fipy.variables.cellVariable import CellVariable
phase = CellVariable(name='phase field', mesh=mesh)
temperature = CellVariable(name='Temperature', mesh=mesh, value=initial_temperature,hasOld=1)
## The orientation is initialized to a uniform value to denote the# randomly oriented liquid phase## .. raw:: latex##    \IndexClass{ModularVariable}#    \IndexModule{numerix}#    \IndexConstant{\pi}{pi}## ..#from fipy.variables.modularVariable import ModularVariable
from fipy.tools.numerix import pi
theta = ModularVariable(
    name='theta',
    mesh=mesh,
    value=-pi + 0.0001,
    hasOld=1
    )## Four different solid circular domains are created at each corner of# the domain with appropriate orientations#L = dx*N
LL = L/60.0
x, y = mesh.getCellCenters()[...,0], mesh.getCellCenters()[...,1]

import random
num_grains = 1
list_of_grains = []
for grain in range(num_grains):
    xx = L*random.random()
    yy = L*random.random()
    angle = 2*pi*random.random() - pi
    list_of_grains.append((xx,yy,angle))
for a, b, thetaValue in list_of_grains:
    segment = (x - a)**2 + (y - b)**2 < (LL)**2
    phase.setValue(1., where=segment)
    theta.setValue(thetaValue, where=segment)
## The `phase` equation is built in the following way. The source term is# linearized in the manner demonstrated in `examples.phase.simple.input`# (Kobayashi, semi-implicit). Here we use a function to build the equation,# so that it can be reused later.## .. raw:: latex##    \IndexClass{TransientTerm}#    \IndexClass{ExplicitDiffusionTerm}#    \IndexClass{ImplicitSourceTerm}## ..#from fipy.terms.transientTerm import TransientTerm
from fipy.terms.explicitDiffusionTerm import ExplicitDiffusionTerm
from fipy.terms.implicitSourceTerm import ImplicitSourceTerm
#def buildPhaseEquation(phase, theta):

    mPhiVar = phase - 0.5 + temperature * phase * (1 - phase)
    thetaMag = theta.getOld().getGrad().getMag()
    implicitSource = mPhiVar * (phase - (mPhiVar < 0))
    implicitSource += (2 * s + epsilon**2 * thetaMag) * thetaMag

    return TransientTerm(phaseTransientCoeff) == \
              ExplicitDiffusionTerm(alpha**2) \
              - ImplicitSourceTerm(implicitSource) \
              + (mPhiVar > 0) * mPhiVar * phase
#phaseEq = buildPhaseEquation(phase, theta)
## The `theta` equation is built in the following way. The details for# this equation are fairly involved, see J.A. Warren *et al.*. The main# detail is that a source must be added to correct for the# discretization of `theta` on the circle.  The source term requires the# evaluation of the face gradient without the modular operators.## .. raw:: latex##    \IndexModule{numerix}#    \IndexFunction{exp}## ..#from fipy.tools.numerix import exp
def buildThetaEquation(phase, theta):

    phaseMod = phase + ( phase < thetaSmallValue ) * thetaSmallValue
    phaseModSq = phaseMod * phaseMod
    expo = epsilon * beta * theta.getGrad().getMag()
    expo = (expo < 100.) * (expo - 100.) + 100.
    pFunc = 1. + exp(-expo) * (mu / epsilon - 1.)

    phaseFace = phase.getArithmeticFaceValue()
    phaseSq = phaseFace * phaseFace
    gradMag = theta.getFaceGrad().getMag()
    eps = 1. / gamma / 10.
    gradMag += (gradMag < eps) * eps
    IGamma = (gradMag > 1. / gamma) * (1 / gradMag - gamma) + gamma
    diffusionCoeff = phaseSq * (s * IGamma + epsilon**2)

    thetaGradDiff = theta.getFaceGrad() - theta.getFaceGradNoMod()
    sourceCoeff = (diffusionCoeff * thetaGradDiff).getDivergence()

    from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
    return TransientTerm(thetaTransientCoeff * phaseModSq * pFunc) == \
               ImplicitDiffusionTerm(diffusionCoeff) \
               + sourceCoeff
## .. raw:: latex##    \IndexClass{ImplicitDiffusionTerm}## ..#thetaEq = buildThetaEquation(phase, theta)
## If the example is run interactively, we create viewers for the phase# and orientation variables. Rather than viewing the raw orientation,# which is not meaningful in the liquid phase, we weight the orientation# by the phase## .. raw:: latex##    \IndexModule{viewers}## ..######## conduction and transient temperature ############################################
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
thermal_conductivity = 80.4  # W/mk
density = 7800 # kg/m^3
heat_capacity = 449.45 # j/(kg*K)
heat_of_fusion = 247292 # j/kg

#FiPy hates division, it burns the skin, so we introduce some variables so it doesnt have to divied
#A=1.0/(density*heat_capacity)
#B=heat_of_fusion/(heat_capacity*timeStepDuration)
A= 2.852e-7
B= 27510.5

temperature_Eq = TransientTerm() == ImplicitDiffusionTerm(thermal_conductivity*A) +  (phase-phase.getOld)*B # <---- this last term causes the problem 





######################################################################
if __name__ == '__main__':
    from fipy.viewers import make
    phaseViewer = make(vars=phase,
                       limits={'datamin': 0., 'datamax': 1.})
    thetaProd = -pi + phase * (theta + pi)
    thetaPhase = phase*theta
    thetaProductViewer = make(vars=thetaProd,
                              limits={'datamin': -pi, 
                                      'datamax': pi})
    phaseViewer.plot()
    thetaProductViewer.plot()
    from fipy.viewers.tsvViewer import TSVViewer
    thetaSave = TSVViewer(vars = thetaProd)
    phaseSave = TSVViewer(vars = phase)
    target = "/Users/edwin/Classes/MSE597I/FiPyExamples/GrowthExample/"
    for i in range(steps):
        theta.updateOld()
        phase.updateOld()
        thetaEq.solve(theta, dt=timeStepDuration)
        phaseEq.solve(phase, dt=timeStepDuration)
        temperature_Eq.solve(temperature, dt=timeStepDuration)
        print i
        if __name__ == '__main__':
            phaseViewer.plot()
            thetaProductViewer.plot()
            save = False
            if save is True:
                if i <10:
                    thetaProductViewer.plot(filename = target +  "theta00%i.eps"%i)                    thetaSave.plot(filename = target + "thetafield00%i.txt"%i)                    phaseSave.plot(filename = target + "phasefield00%i.txt"%i)                elif i<100:                    thetaProductViewer.plot(filename = target + "theta0%i.eps"%i)                    thetaSave.plot(filename = target + "thetafield0%i.txt"%i)                    phaseSave.plot(filename = target + "phasefield0%i.txt"%i)                else:                    thetaProductViewer.plot(filename = target + "theta%i.eps"%i)                    thetaSave.plot(filename = target + "thetafield%i.txt"%i)                    phaseSave.plot(filename = target + "phasefield%i.txt"%i)