Re: Seeing solution on variable mesh

[Daniel Wheeler]

Hi Terry, Answers below. On Aug 5, 2006, at 4:36 PM, [EMAIL PROTECTED] wrote: Dan,     I've been poking around with my 1D thermal problem.  The main change that I am trying to make is to expand the mesh GREATLY around the solution region.  You may recall that originally the total span of x was about 270 nm.   Now, I want to add material on the low and high x sides so that the original thin film structure is embedded in the center (roughly) of a 10 mm long model. I have solved problems with a great deal of different mesh sizes, the key is to taper the mesh from course to fine in a sensible way. Basically, adjacent elements should not have massive size differences, it effects accuracy. Also the matrix should not have rows that are numerically insignificant, this can cause problems. I have worked with meshes that vary from 5nm to 1 micron. You will probably be alright with the dimensions you are suggesting since you are in 1D.     This code seems to run without errors.    (See attached file: Yukiko_DW3d.py)      However, the output plot is uninteresting, showing T=0 for all visible x.  I realize that my "interesting" region is initially TINY in the center of the 10 mm wide model, BUT the heat should diffuse outward to a noticeable degree during the solution.  And what I really am after is a reliable peak temperature rise, given this extensive buffer material surrounding the heated zone.  Before, my peak temperature rise was NOT reliable because my T=0 boundaries were far too close to the heated zone.     Can you see why I am not observing any non-zero temperature rise?   Is the mesh size too variable for the solver to get a sensible solution?  Is the plot resolution insufficient to discern anything? Changed a few things and it seems to be doing something. Don't know if it makes sense:

# twm 8/2,5u/06 - change spatial scale to ~10mm
shift = 5.e-3
bot = 5.e-3 + 1.e-6
dx = (1.e-10,1.e-3,1.e-3,1.e-3,1.e-3,.5e-3,.3e-3,.1e-3,.05e-3,.05e-3,3.e-9,3.e-9,3.e-9,3.e-9,3.e-9,3.e-9,3.e-9,20.e-9,60.e-9,200.e-9,500.e-9,3000.e-9,1.e-5,3.e-5,1.e-4,3.e-4,5.e-4,1.e-3,1.e-3,1.e-3,1.e-3,1.e-10)
dt = 2.e-5

ab = 0.412
alpha = 1 / 30.e-9
Fl = 30.e4                
tf = 25.e-3
go = 2. / tf
gxo = ab * alpha * Fl

specificHeat = (0.9268E3, 2.52E6, 3.56E6, 1.88E6, 0.9268E3)
conductivity = (0.02, 2.0, 20.0, 0.7, 0.02)

from fipy.meshes.grid1D import Grid1D
mesh = Grid1D(dx=dx)

from fipy.variables.cellVariable import CellVariable
from fipy.variables.variable import Variable

var = CellVariable(mesh=mesh)
transientCoeff = CellVariable(mesh=mesh)
diffusionCoeff = CellVariable(mesh=mesh)
time = Variable(0)

for i in range(len(specificHeat)):
     mask = mesh.getCellCenters()[...,0]
     transientCoeff.setValue(specificHeat[i], where=mask)
     diffusionCoeff.setValue(conductivity[i], where=mask)

from fipy.terms.transientTerm import TransientTerm
from fipy.terms.diffusionTerm import DiffusionTerm
from fipy.tools import numerix

X = mesh.getCellCenters()[...,0]

Xmask = (X >= shift) * (X <= (shift + bot)) * (time < tf)

source = Xmask * go * gxo * (1 - time / tf) * numerix.exp(-alpha * (X - shift) * Xmask)

##if (X < shift): source = 0. 
##elif ((X >= shift) and (X <= (shift + bot))):
##     source = go * gxo * (X > shift) * (time < tf) * (1 - time / tf) * numerix.exp(-alpha * (X - shift))
##elif (X > (shift + bot)):      
##     source = 0. 

eqn = TransientTerm(transientCoeff) == DiffusionTerm(diffusionCoeff) + source

from fipy.viewers import make
viewer = make(var)
# viewer = make(var, limits={'datamax':100})

from fipy.boundaryConditions.fixedValue import FixedValue
BCs = (FixedValue(faces=mesh.getFacesLeft(), value=0),
       FixedValue(faces=mesh.getFacesRight(), value=0))

while (time < 2.0 * tf).getValue():
     eqn.solve(var, boundaryConditions=BCs, dt=dt)
     print 'X - shift',X - shift
     print 'source',source
#     raw_input("press key")
     time.setValue(time + dt)

     viewer.plot()

raw_input("finished")

You were setting the source term to zero everywhere. The source in this case is an implicitly defined variable that is automatically updated if anything it depends on is changed, "time" in this case. The element wise conditional operations have to be part of the _expression_ in order to work correctly. Also asking if "X > 0" returns an array if X is an array and thus you can't very well use conditionals with an array of ones and zeros.      One tricky thing that I had to do was define the non-zero SOURCE term only in the thin film region, lest I get some numerical overflow (probably from the exp function, I imagine).  I got past that error using your old SHIFT variable, along with my new BOT variable.  But now I don't know if my source is actually non-zero anywhere, given the large dynamic range of my spatial variable x. All you have to do is print its value to see.     >>> print ' my source term',source     Another thing here that is bit awkward for me is the specification of the mesh INTERVALS dx instead of the mesh nodes x.  Can I just as well define the mesh nodes along x directly? Yes, define X, the array of nodal positions, then do,       >>> dx = X[1:] - X[:-1] to define dx. I realize that FiPy wants to use cell CENTERS in the solution, and that's OK.  I would just like to specify my mesh nodes [x1, x2, ...., xn] directly for input and for plotting output. You can ask the mesh for its vertex coordinates;     >>> mesh.getVertexCoords()     Thanks. Terry McDaniel Seagate Research     (209) 295-6735 <Yukiko_DW3d.py> Daniel Wheeler