Daniel Wheeler wrote:
> Hi Michael,
>
> Reply is below.
>
> On Fri, Feb 22, 2008 at 2:13 PM, Michael Waters <[EMAIL PROTECTED]> wrote:
>
>> I am trying to compute an interface angle in the same way that was done
>> in the paper "Extending phase field models of solidification to
>> polycrystalline materials" by James Warren, Ryo Kobayashi, Alexander
>> Lobkovsky, Craig Carter. I used a snippet of code from the fipy
>> anisotropy example (section 9.4), but this returns an error that I'm not
>> sure how to correct. If anyone knows a better way to do this or how to
>> work around the error, I'd love to hear.
>>
>> This is the code that I'm using, I'd like it if I could keep theta as a
>> modular variable:
>>
>> phase = CellVariable(name='phase field', mesh=mesh, value = 0.0)
>> theta = ModularVariable( name='theta', mesh=mesh, value=-pi + 0.0001,
>> hasOld=1 )
>>
>> phaseY = phase.getFaceGrad().dot((0, 1))*1.0
>> phaseX = phase.getFaceGrad().dot((1, 0))*1.0
>> psi= arctan2(phaseY, phaseX) + theta
>>
>
> I think this error is simply due to doing mixed operations with face
> and cell values. It has
> nothing to do with the fact that theta is a ModularVariable. The error
> message is not particularly helpful though.
>
> To fix it do the following. Make phaseY and phaseX cell values
>
> phaseY = phase.getGrad().dot((0, 1))
> phaseX = phase.getGrad().dot((1, 0))
>
> Then
>
> psi= arctan2(phaseY, phaseX) + theta
>
> works. FYI is not a ModularVariable. Do you want it to be?
>
> Looking into this I realize that ModularVariable does not have
> addition defined. This doesn't work
>
> >>> print theta + theta
> [-6.28298531 -6.28298531 -6.28298531 -6.28298531] 1
>
> This does work
>
> >>> print theta - (-theta)
> [ 0.0002 0.0002 0.0002 0.0002] 1
>
> which I believe is correct. Also this works
>
> >>> print 4 * theta - theta
> [-3.14129265 -3.14129265 -3.14129265 -3.14129265] 1
>
> I am not sure why the addition operator was never included. Probably
> there is a reason, but I can't remember.
> I am now worried that the other operators need to be fixed. What about this?
>
> >>> print theta * 4
> [-12.56597061 -12.56597061 -12.56597061 -12.56597061] 1
>
> Anyhow, if you are only doing limited operations, then I think you can
> work around it.
>
> Is ModularVariable worth fixing up to work properly? Or is there a
> better way such as just defining a calculation
> and then just using a mod operator like this
>
> >>> a = CellVariable(..)
> >>> b = CellVariable(..)
> >>> c = (a + b) % numerix.pi
>
> Seems more natural to me. The % operator should return a BinOp, which
> is currently not the case (Is it?).
> We would then have to do something about the weird face grads that
> need to be done in some of
> the examples in a modular fashion, but they probably could be defined
> as special variables or even in the script.
> I think this makes much more sense. Of course, the following,
>
> >>> d = ((a + b) % numerix.pi) + a
>
> would not be a modular variable, which also makes sense to me.
>
> Cheers
>
> --
> Daniel Wheeler
>
> >Thank you, Daniel Wheeler.
I'm running into a new error however, and was hoping someone could
explain to me what is happening. I am trying to create a simple
solidification model.
The problem comes up as this:
TypeError: unsupported operand type(s) for +: 'binOp' and 'unOp'
I am new to fipy and python so im having trouble conceptualizing binOp
vs. unOp. The error comes up around line 120.
thanks!
-Mike W.
#!/usr/bin/env python
from fipy.tools.parser import parse
numberOfElements = parse('--numberOfElements', action = 'store',
type = 'int', default = 40)
numberOfSteps = parse('--numberOfSteps', action = 'store',
type = 'int', default = 1000)
steps = numberOfSteps
from fipy.tools.numerix import sqrt
N = int(sqrt(numberOfElements))
L = 2.5 * N / 100.
dL = L / N
dx = L/N
dy = L/N
from fipy.meshes.grid2D import Grid2D
mesh = Grid2D(dx=dx, dy=dx, nx=N, ny=N)
timeStepDuration = 1e-5
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= -1.0
Boundary_Temperature = -1.0
#############################
tau = 3e-4
c = 0.02
Fold_Symetry = 4.0
kappa1 = 0.9
kappa2 = 20.0
tempDiffusionCoeff = 2.25
from fipy.variables.cellVariable import CellVariable
phase = CellVariable(name='phase field', mesh=mesh, hasOld = 1)
temperature = CellVariable(name='Temperature', mesh=mesh,
value=initial_temperature,hasOld=1)
from fipy.variables.modularVariable import ModularVariable
from fipy.tools.numerix import pi
theta = ModularVariable(
name='theta',
mesh=mesh,
value=-pi + 0.0001,
hasOld=1
)
#
#
L = dx*N
LL = dx*5
x, y = mesh.getCellCenters()[...,0], mesh.getCellCenters()[...,1]
import random
num_grains = 5
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)
#
# .. raw:: latex
from fipy.tools.numerix import pi, arctan, arctan2, tan
phaseY = phase.getGrad().dot((0, 1))
phaseX = phase.getGrad().dot((1, 0))
psi= arctan2(phaseY, phaseX) + theta
Phi = tan(Fold_Symetry * psi / 2.0)
PhiSq = Phi**2
beta_2 = (1.0 - PhiSq) / (1.0 + PhiSq)
beta_2Psi = -Fold_Symetry * 2 * Phi / (1.0 + PhiSq)
A = alpha**2 * c * (1.+ c * beta_2) * beta_2Psi
D = alpha**2 * (1.+ c * beta)**2
#############
dxi = (A*phase).getGrad()._take((1, 0), axis = 1) * (-1, 1)
anisotropySource =(dxi).getDivergence()
###########
########
#
print 'building phase equation'
from fipy.terms.transientTerm import TransientTerm
from fipy.terms.explicitDiffusionTerm import ExplicitDiffusionTerm
from fipy.terms.implicitSourceTerm import ImplicitSourceTerm
#
def buildPhaseEquation(phase, theta):
mVar = phase - 0.5 - kappa1 / pi * arctan(kappa2 * temperature)
thetaMag = theta.getOld().getGrad().getMag()
BB= mVar * ((mVar < 0.0) - phase) + (mVar > 0.0) * mVar * phase #variable
BB is really meant to reduce the line length below
print 'BB + anisotropySource =', thetaMag + anisotropySource ### if this
line is enabled, the error returned is :
##### TypeError: unsupported operand type(s) for +: 'unOp' and 'unOp'
## this is odd because the addition below causes this error:
#### TypeError: unsupported operand type(s) for +: 'binOp' and 'unOp'
#to me, it seems that these should be the same type just in different
lines.
return TransientTerm(phaseTransientCoeff) == ImplicitSourceTerm(BB+
anisotropySource) + ExplicitDiffusionTerm(D) ### anisotropySource term is
causing the problem.
phaseEq = buildPhaseEquation(phase, theta)
#
print 'theta equation building'
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()
print phaseFace
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
thetaEq = buildThetaEquation(phase, theta)
#
from fipy.solvers.linearLUSolver import LinearLUSolver
solver = LinearLUSolver(tolerance=1e-10)
#
####### conduction and transient temperature
############################################
print 'building thermal diffusion equation'
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
tempDiffusionCoeff = 2.25
temperature_Eq = TransientTerm() == ImplicitDiffusionTerm(tempDiffusionCoeff) +
(phase-phase.getOld())/timeStepDuration
from fipy.boundaryConditions.fixedValue import FixedValue
BCS_temperature=(FixedValue(mesh.getFacesRight(),Boundary_Temperature),FixedValue(mesh.getFacesLeft(),Boundary_Temperature),FixedValue(mesh.getFacesTop(),Boundary_Temperature),FixedValue(mesh.getFacesBottom(),Boundary_Temperature))
######################################################################
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})
temp_view = make(vars=temperature,
limits={'datamin': Boundary_Temperature,
'datamax': 0.0 })
phaseViewer.plot()
thetaProductViewer.plot()
temp_view.plot()
from fipy.viewers.tsvViewer import TSVViewer
thetaSave = TSVViewer(vars = thetaProd)
phaseSave = TSVViewer(vars = phase)
Temperature_Save = TSVViewer(vars = temperature)
target = "/Users/mjwaters/DVLP/Solidify"
print 'running model'
for i in range(steps):
theta.updateOld()
phase.updateOld()
temperature.updateOld()
thetaRES=1.0
phaseRES=1.0
temperatureRES=1.0
while temperatureRES > 1e-3 or phaseRES > 1e-3 or thetaRES > 1e-3:
thetaRES= thetaEq.sweep(theta,solver = solver ,
dt=timeStepDuration)
phaseRES= phaseEq.sweep(phase,solver = solver ,dt=timeStepDuration)
temperatureRES=temperature_Eq.sweep(temperature,solver = solver,
dt=timeStepDuration, boundaryConditions = BCS_temperature)
print thetaRES, phaseRES, temperatureRES
print i
print 'temperature max = ', max(temperature)
if __name__ == '__main__':
phaseViewer.plot()
thetaProductViewer.plot()
temp_view.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)
