Hi Drew, The end problem that you'd like to solve is quite complicated and there is a lot to consider. My recommendation would be to not use FiPy to get started, but just use Scikit-fmm, http://pythonhosted.org/scikit-fmm/, to understand the level set method and velocity extensions as well as advecting the interface. FiPy uses Scikit-fmm for its level set capability. You can do this in 1D explicitly and implement the bulk equations using a 1D finite difference with Numpy. This well help you understand the numerical implementation of the interface velocity from the bulk equations, which can be quite difficult. FiPy does not take care of a general interfacial boundary condition so you really need to understand how to implement this numerically. Once you have something you understand and can compare with, then at that point it might be a good idea to start using FiPy for a fully implicit implementation of the bulk equations.
Hope that helps, Daniel On Mon, Feb 13, 2017 at 12:33 PM, Drew Davidson <davidson...@gmail.com> wrote: > Hi, > > > My ultimate goal: a numerical model very similar to doi:10.2298/SOS1001033N > (Computer simulation of rapid solidification with undercooling: A case study > of spherical ceramics sample on metallic substrate). The full text pdf of > that paper is easily found online. I'm interested in solidification of a > pure metal with interface velocity a function of interface temperature > (thermally undercooled solidification of a pure metal). My knowledge of > solidification, Python, and FiPy is quite basic. I just discovered FiPy a > few weeks ago. > > > I want to start by considering the level set method in FiPy. Start with the > basic 1D Stefan problem: solidification of a pure material on a finite > domain. I know nothing about the phase field method and whether it is > applicable to the ultimate goal in paragraph 1, in which interface velocity > is a specified function of interface temperature. > > > I found a thread from 2007 > https://www.mail-archive.com/fipy@nist.gov/msg00320.html, but the code of > Daniel Wheeler is too advanced for me to understand, and does not run > unmodified in current FiPy. Without being able to comprehend this code in > 1D, extension to 2D cylindrical coordinates seems unattainable. > > > I looked at the level set example of a 'simple trench system:' > http://www.ctcms.nist.gov/fipy/examples/levelSet/generated/examples.levelSet.electroChem.howToWriteAScript.html > > > I tried to imitate this example (the metal ion Cm part) with the scripts > below, but am stuck as to how to couple the interface velocity into my > problem as an unknown. In my script below, the interface velocity appears to > remain at the initial specified value, since it apparently never gets > updated or coupled into the system. It would seem necessary to have the > Stefan condition as the coupling equation, but I am unable to look at the > FiPy manual and FiPy Python code and identify the necessary commands. I > looked for the equivalent of the Stefan condition (not expressed as a source > term, but expressed in terms of gradients of metal ion concentration on > either side of the interface) in the code for the example Simple Trench > System, but was unable to identify it. The Stefan condition is currently > represented in my script as a source term, but I do not see how to couple > Vinterface into the unknowns. It seems like one approach could be to replace > Vinterface in the source term with an expression involving gradients of > temperature on either side of the interface, but I am stuck as to how to > write this in the language of FiPy in both 1D and eventually 2D (first > Cartesian 2D then cylindrical 2D). > > > I would love to see code that is can be obviously generalized to 2D > Cartesian then 2D Cylindrical, instead of being specialized only to 1D. The > code in https://www.mail-archive.com/fipy@nist.gov/msg00320.html might be > 1D-only and/or non-obvious to translate to 2D, for example. > > > Thanks > > > MY MAIN PYTHON SCRIPT > > > > # from IPython import embed #then put embed() in your script at a point > where you want to be at ipython prompt > > > # take the approach of having variables with conventional names e.g. > temperature, then set to 1.0 etc for nondimensional > > # this lets you choose different nondimensionalization schemes > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > nx = 20 > > L=1. > > cellSize = L/nx > > #there is no use for a graded mesh since grid is fixed and interface is > moving > > mesh = Grid1D(nx=nx, Lx=L) > > x=mesh.cellCenters > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > numberOfCellsInNarrowBand = 10 > > narrowBandWidth = numberOfCellsInNarrowBand * cellSize > > cflNumber=.2 > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > #level set method distance variable, initialized to -1. > > distanceVar = DistanceVariable(name = 'distance variable',mesh = mesh,value > = -1.,hasOld = 1) > > #TODO I think the interface has to advance a little into the domain to > establish where it is? > > #alternative is to have extra cells to the left of the actual x=0, as in > simpleTrenchSystem.py > > distanceVar.setValue(1., where=(x > cellSize)) > > distanceVar.calcDistanceFunction(order=2) > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > #thermal diffusivity > > alphaSolid=1. > > alphaLiquid=1. > > > # #thermal conductivity > > # kSolid=1. > > # kLiquid=1. > > > # #mass density > > # rhoSolid=1. > > # rhoLiquid=1. > > > #specific heat > > cpSolid=1. > > cpLiquid=1. > > > latentHeat=1. > > Tm=.5 #melting T > > > temperature = > CellVariable(name="temperature",mesh=mesh,value=1.,hasOld=False) > > > #left boundary is step change in temperature > > valueLeft=0. > > temperature.constrain(valueLeft, mesh.facesLeft) > > > #right boundary is zero flux, which is default BC in FiPy > > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > #TODO this scheme is not going to work. interface velocity starts out with > initial guess. how is it updated or linked to anything? > > #simple trench electro problem: interface velocity was a function of metal > ion bulk concentration and theta > > > interfaceVelocVar=CellVariable(name='interface > velocity',mesh=mesh,value=0.1) #TODO initializing V at zero results in > errors > > > advectionEquation = TransientTerm() + AdvectionTerm(interfaceVelocVar) > > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > energyEquation=buildEnergyEquation(temperature=temperature,distanceVar=distanceVar,interfaceVelocVar=interfaceVelocVar,cellSize=cellSize,alphaSolid=alphaSolid,alphaLiquid=alphaLiquid,latentHeat=latentHeat,cpSolid=cpSolid,cpLiquid=cpLiquid,Tm=Tm) > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > > viewer = Viewer(vars=(temperature,),datamin=0., datamax=1.) > > viewer.plot() > > raw_input() > > > levelSetUpdateFrequency = int(0.8 * narrowBandWidth / (cellSize * cflNumber > * 2)) #TODO where does this come from? > > > numberOfSteps=80 > > > for step in range(numberOfSteps): > > if step>5 and step % 5 == 0 and viewer is not None: > > viewer.plot() > > raw_input() > > > if step % levelSetUpdateFrequency == 0: > > distanceVar.calcDistanceFunction(order=2) > > > distanceVar.updateOld() > > > distanceVar.extendVariable(interfaceVelocVar, order=2) > > dt = cflNumber * cellSize / interfaceVelocVar.max() > > > advectionEquation.solve(distanceVar, dt = dt) > > energyEquation.solve(temperature,dt=dt) > > > > > > > > MY SUBROUTINE SCRIPT TO BUILD ENERGY EQUATION > > > from fipy.terms.implicitSourceTerm import ImplicitSourceTerm > > from fipy.terms.transientTerm import TransientTerm > > from fipy import FaceVariable, CellVariable, DiffusionTerm > > > def > buildEnergyEquation(temperature,distanceVar,interfaceVelocVar,cellSize,alphaSolid=1.0,alphaLiquid=1.0,latentHeat=1.0,cpSolid=1.,cpLiquid=1.,Tm=0.): > > r''' > > #based on metalIonDiffusionEquation.py from NIST > > > #solve the one dimensional solidification problem from materials processing > class, > > #9/28/1999 lecture: solidification in finite domain; page 105 in class notes > > #level set method > > ''' > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > mesh = distanceVar.mesh > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > #http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html > > #trying to adapt the position-dependent diffusivity example to > distanceVar-dependent diffusivity > > D = FaceVariable(mesh=mesh, value=alphaSolid) > > > #https://www.mail-archive.com/fipy@nist.gov/msg00854.html > > distVarAtFaceCenters=distanceVar.getHarmonicFaceValue() > > > D.setValue(alphaLiquid, where=(distVarAtFaceCenters > 0.)) > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > > cp=CellVariable(name='specific heat',mesh=mesh,value=cpSolid,hasOld=0) > > cp.setValue(cpLiquid, where=(distanceVar > 0)) > > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > > #TODO is it OK that interfaceVelocVar is a CellVariable? in > metalIonDiffusion.py, the default interface V was a scalar > > #here is where the Stefan condition is expressed as a source term > > #dividing by solution variable because of ImplicitSourceTerm: see FiPy > manual > > coeff = interfaceVelocVar * latentHeat * distanceVar.cellInterfaceAreas / > (mesh.cellVolumes) / cp / temperature > > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > #need to fix interface T at solidification T > > # > http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-internal-boundary-conditions > > mask1=((distanceVar >= -cellSize/2.) & (distanceVar<cellSize/2.)) #TODO will > this always select a single cell center? is syntax correct? ipython > > largeValue=1e+10 > > > # > (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( > > #TODO or the == version? > > #TODO check signs > > transientCoeff=1. > > eq = TransientTerm(transientCoeff) - DiffusionTerm(D) - > ImplicitSourceTerm(coeff) + ImplicitSourceTerm(largeValue*mask1) - > largeValue*mask1*Tm > > > return eq > > > > > END OF MESSAGE > > > _______________________________________________ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > -- Daniel Wheeler _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
