Dear FiPy group,
I am trying to solve a very simple Cahn-Hilliard equation in FiPy, where I
am using a scaling parameter, L, to adjust the dimensions of the domain.
When I set the length, L = 40, the equation cannot be solved (test1.py).
However, t if I set L = 0.4 e-5.*Magnification
where
Magnification = 1.0e7
which numerically is ALSO equal to 40, the equation becomes
solvable(test2.py). I am baffled by this, and cannot understand why one
simulation works and the other one doesn't.
By typing,
diff test1.py test2.py
in a shell terminal you can find:
16c16,17
< L = 40 # dimension of simulation domain
---
> Magnification = 1.0e7 # magnification
> L = 0.4e-5*Magnification # dimension of simulation domain
which proves that the both have the same L value. The only difference is
how it is being inputed. I am attaching both scripts, test1.py and test2.py
I would appreciate your help.
Thanks,
Rui
########### IMPORT STATEMENTS #############
from fipy import CellVariable, DistanceVariable, buildAdvectionEquation, Viewer, Grid2D, dump
import numpy as np
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.transientTerm import TransientTerm
from fipy.tools.numerix import dot
from fipy.solvers.pysparse.linearPCGSolver import LinearPCGSolver
import matplotlib.pyplot as plt
from matplotlib import cm, pylab
import time
from fipy.boundaryConditions.fixedFlux import FixedFlux
from fipy.boundaryConditions.nthOrderBoundaryCondition import NthOrderBoundaryCondition
N = 40 # mesh size 40X40
nsteps = 15000 # total number of time steps to run
L = 40 # dimension of simulation domain
dt = 0.025 # time step [s]
R = 8.314 # gas constant [J/K/mol]
Temperature = 1200.0 # chamber temperature [K]
m = Grid2D(nx=N, ny=N, Lx=L, Ly=L) # define the mesh
######### Define and initialize concentration variables ########
concentration = CellVariable(mesh=m, name='concentration') # define concentration variable, marks the solute concentration.
def noise(length,percentage):
noise = 1 + (np.random.random(length)-0.5)*2*percentage
return noise
concentration.setValue(0.5)
concentration[:]=concentration.value*noise(N**2, 0.30)
############# build the diffusion equation inside the solid phase ##########
solidPhaseDerivative = 0.5*(12.0*concentration**2 - 12.0*concentration + 2.0) # the second derivative of the energy well
kappa = 1.0 # energy coefficient along concentration gradient
M_solid = 1.0 # solid phase diffusivity [m^2/s]
diffTerm2 = ImplicitDiffusionTerm(coeff = (M_solid*solidPhaseDerivative,))
diffTerm4 = ImplicitDiffusionTerm(coeff = (M_solid, kappa))
eqch = TransientTerm() - diffTerm2 + diffTerm4
solver = LinearPCGSolver(tolerance=1.0e-10, iterations=1000)
##### define the boundary condition###
BCs = (FixedFlux(m.facesRight, 0),
FixedFlux(m.facesLeft, 0),
NthOrderBoundaryCondition(m.facesLeft, 0, 3),
NthOrderBoundaryCondition(m.facesRight, 0, 3),
NthOrderBoundaryCondition(m.facesTop, 0, 3),
NthOrderBoundaryCondition(m.facesBottom, 0, 3))
for step in range(nsteps):
################ plot the cocentration variable every 50 steps ##########
if step%50==0:
values = (concentration.value).reshape((N,N))
fig = plt.figure()
fig.frameon=False
ax = fig.add_subplot(111);im=plt.imshow(values,vmin=0.0,vmax=1.0,cmap=cm.jet)
ax1=pylab.gca();plt.subplots_adjust(left=0,right=1.0,top=1.0,bottom=0.)
ax1.linewidth=4.;pylab.gca().set_xticks(());pylab.gca().set_yticks(())
pylab.gca().spines['left'].set_linewidth(0.);pylab.gca().spines['right'].set_linewidth(0.);
pylab.gca().spines['top'].set_linewidth(0.);pylab.gca().spines['bottom'].set_linewidth(0.);
im.set_interpolation('bicubic')
plt.savefig('output1/plot_%.5d.pdf' % step)
####################solve equations#############################
eqch.solve(concentration, boundaryConditions=BCs, solver=solver, dt=dt)
########### IMPORT STATEMENTS #############
from fipy import CellVariable, DistanceVariable, buildAdvectionEquation, Viewer, Grid2D, dump
import numpy as np
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.transientTerm import TransientTerm
from fipy.tools.numerix import dot
from fipy.solvers.pysparse.linearPCGSolver import LinearPCGSolver
import matplotlib.pyplot as plt
from matplotlib import cm, pylab
import time
from fipy.boundaryConditions.fixedFlux import FixedFlux
from fipy.boundaryConditions.nthOrderBoundaryCondition import NthOrderBoundaryCondition
N = 40 # mesh size 40X40
nsteps = 15000 # total number of time steps to run
Magnification = 1.0e7 # magnification
L = 0.4e-5*Magnification # dimension of simulation domain
dt = 0.025 # time step [s]
R = 8.314 # gas constant [J/K/mol]
Temperature = 1200.0 # chamber temperature [K]
m = Grid2D(nx=N, ny=N, Lx=L, Ly=L) # define the mesh
######### Define and initialize concentration variables ########
concentration = CellVariable(mesh=m, name='concentration') # define concentration variable, marks the solute concentration.
def noise(length,percentage):
noise = 1 + (np.random.random(length)-0.5)*2*percentage
return noise
concentration.setValue(0.5)
concentration[:]=concentration.value*noise(N**2, 0.30)
############# build the diffusion equation inside the solid phase ##########
solidPhaseDerivative = 0.5*(12.0*concentration**2 - 12.0*concentration + 2.0) # the second derivative of the energy well
kappa = 1.0 # energy coefficient along concentration gradient
M_solid = 1.0 # solid phase diffusivity [m^2/s]
diffTerm2 = ImplicitDiffusionTerm(coeff = (M_solid*solidPhaseDerivative,))
diffTerm4 = ImplicitDiffusionTerm(coeff = (M_solid, kappa))
eqch = TransientTerm() - diffTerm2 + diffTerm4
solver = LinearPCGSolver(tolerance=1.0e-10, iterations=1000)
##### define the boundary condition###
BCs = (FixedFlux(m.facesRight, 0),
FixedFlux(m.facesLeft, 0),
NthOrderBoundaryCondition(m.facesLeft, 0, 3),
NthOrderBoundaryCondition(m.facesRight, 0, 3),
NthOrderBoundaryCondition(m.facesTop, 0, 3),
NthOrderBoundaryCondition(m.facesBottom, 0, 3))
for step in range(nsteps):
################ plot the cocentration variable every 50 steps ##########
if step%50==0:
values = (concentration.value).reshape((N,N))
fig = plt.figure()
fig.frameon=False
ax = fig.add_subplot(111);im=plt.imshow(values,vmin=0.0,vmax=1.0,cmap=cm.jet)
ax1=pylab.gca();plt.subplots_adjust(left=0,right=1.0,top=1.0,bottom=0.)
ax1.linewidth=4.;pylab.gca().set_xticks(());pylab.gca().set_yticks(())
pylab.gca().spines['left'].set_linewidth(0.);pylab.gca().spines['right'].set_linewidth(0.);
pylab.gca().spines['top'].set_linewidth(0.);pylab.gca().spines['bottom'].set_linewidth(0.);
im.set_interpolation('bicubic')
plt.savefig('output2/plot_%.5d.pdf' % step)
####################solve equations#############################
eqch.solve(concentration, boundaryConditions=BCs, solver=solver, dt=dt)
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
