is it easy to get fipy to make full use of dual CPUs in a single workstation?
Hello, Does anybody know how easy it is to get fipy installed such that it makes use of both CPUs in a dual CPU workstation? The OS would be Ubuntu or OpenSUSE. Will fipy use both CPUs out of the box? Or would I need some sophistication to get it to work? I know nothing about Trilinos etc. I don't currently have the workstation. I was somewhat interested in buying something like the used workstation on the ithaca ny craigslist (specs given below) since it seems like it might be able to do fipy phase field calculations for larger problems because of what seems like relatively numerous cores and threads, but I would not want to buy a such a device, if there was a chance I would be unable to configure fipy to make full use of it. I have already observed that with my single CPU system (i5-3550p), fipy seems to max out its four cores in ubuntu and opensuse without any additional setup. Specs for the Ithaca Craigslist example workstation: Dual Xeon E5-2680 2.7 GHz / 3.5 GHz Turbo (8 cores, 16 with HT each or 16 cores, 32 with HT total) 64GB (8 x 8GB) PC3 Registered ECC Memory This is something I have been curious about for a long time... Thanks. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to apply convection boundary condition in 2D?
Hi Dr. Guyer, I updated my code with your changes. It now matches the analytical solution pretty well, after several grid refinements. This is of course not a comprehensive test. The Python script: https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/ConvectionTestProblem2D_01_20181215Version.py A report showing agreement between FiPy and analytical solution: https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/Report20181216.pdf The README.md of my repo has been updated to provide a better explanation of what is going on. I apologize for not learning git/github well enough to have github properly track changes between us as we went along. Instead, I put comments at the top of my scripts to document how I was looking at your code and documentation. Thanks On Fri, Dec 14, 2018 at 8:20 PM Guyer, Jonathan E. Dr. (Fed) via fipy < fipy@nist.gov> wrote: > Drew - > > I had found the same sign error in the denominator. You're also correct > that the units didn't agree. Multiplying by the face areas is redundant to > taking the divergence. > > I've updated our notes: > > > https://github.com/usnistgov/fipy/blob/dd3420fb71884d74850051ad2280bff525301824/documentation/USAGE.rst > > and your latest code: > > > https://github.com/guyer/Convection2DFiPyExample01/commit/ca59e8955319545f9965705c8fadfcbb5abd5951 > > With these changes, the results look much more plausible. > > - Jon > > > On Dec 12, 2018, at 1:22 PM, Drew Davidson > wrote: > > > > Hi Dr. Guyer, > > > > I looked over your more recent documentation: > > > https://github.com/usnistgov/fipy/blob/70b72b7abb267c85ada47886b8b3573e1819fffc/documentation/USAGE.rst > > > > I am embarrassed I did not understand how to choose values for a and b > in the Robin condition before. > > > > I also cloned your repo to my local computer and ran it: > > https://github.com/guyer/Convection2DFiPyExample01 > > > > When I run the code in your repo, the viewer does show a variation in T > with respect to r, but the size of that variation is extremely tiny. > Essentially, the solution is pinned at the initial temperature. The > variation also has the wrong sign (solid object is warming up), as the > ambient air is colder than the solid object (T_infinity < T_initial), and > the solid object should be cooling down. Maybe this is due to confusion I > sowed by having var be T minus T_infinity; sorry. > > > > I tried another script in my repo: > > > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/ConvectionTestProblem2D_01_20181211Version.py > > > > This script is upgraded in that now, the log file contains information > about the solution (max and min temperatures). I obtain max and min > temperatures using var.value.min() and var.value.max() in order to rule out > goof-ups in my get_solution_along_ray function. > > > > The solution variable var is now exactly the temperature, in order to > get rid of the previous source of confusion in which var was T minus > T_infinity. > > > > The solution in my script is still pinned at the initial temperature. > This persists even if I crank up the convection coefficient to a huge value > (keep multiplying it by 1000 over and over and still seeing solution pinned > at initial condition). > > > > I made handwritten notes on top of a pdf I made via Kile Editor of your > more recent documentation: > > > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsage20181211_Annot_CTTC.pdf > > > > It seems to me that the units of the source terms should match the units > of the terms in the partial differential equation for heat conduction. I > thought that is how it worked from the FiPy examples, but I haven’t gone > back and made sure. I find that for the first source term RobinCoeff*g, the > units don’t seem to match the units of the terms in the PDE for heat > conduction (I didn’t check the 2nd source term since it should be the same > situation). I did this in kind of a hurry so maybe I am goofing up. I have > not studied the finite volume method, and am just going on a hunch. > > > > It also seemed to me that there might be a minus sign missing in the > RobinCoeff. The solution I get is still pinned at initial condition > regardless of a minus sign or no minus sign there. > > > > Thanks > > > > On Fri, Dec 7, 2018 at 6:05 PM Guyer, Jonathan E. Dr. (Fed) via fipy < > fipy@nist.gov> wrote: > > Drew - > > > > > > Thanks to your notes, I found a couple of errors in our Robin > discussion. As you note
Re: how to apply convection boundary condition in 2D?
Hi Dr. Guyer,
I looked over your more recent documentation:
https://github.com/usnistgov/fipy/blob/70b72b7abb267c85ada47886b8b3573e1819fffc/documentation/USAGE.rst
I am embarrassed I did not understand how to choose values for a and b in
the Robin condition before.
I also cloned your repo to my local computer and ran it:
https://github.com/guyer/Convection2DFiPyExample01
When I run the code in your repo, the viewer does show a variation in T
with respect to r, but the size of that variation is extremely tiny.
Essentially, the solution is pinned at the initial temperature. The
variation also has the wrong sign (solid object is warming up), as the
ambient air is colder than the solid object (T_infinity < T_initial), and
the solid object should be cooling down. Maybe this is due to confusion I
sowed by having var be T minus T_infinity; sorry.
I tried another script in my repo:
https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/ConvectionTestProblem2D_01_20181211Version.py
This script is upgraded in that now, the log file contains information
about the solution (max and min temperatures). I obtain max and min
temperatures using var.value.min() and var.value.max() in order to rule out
goof-ups in my get_solution_along_ray function.
The solution variable var is now exactly the temperature, in order to get
rid of the previous source of confusion in which var was T minus T_infinity.
The solution in my script is still pinned at the initial temperature. This
persists even if I crank up the convection coefficient to a huge value
(keep multiplying it by 1000 over and over and still seeing solution pinned
at initial condition).
I made handwritten notes on top of a pdf I made via Kile Editor of your
more recent documentation:
https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsage20181211_Annot_CTTC.pdf
It seems to me that the units of the source terms should match the units of
the terms in the partial differential equation for heat conduction. I
thought that is how it worked from the FiPy examples, but I haven’t gone
back and made sure. I find that for the first source term RobinCoeff*g, the
units don’t seem to match the units of the terms in the PDE for heat
conduction (I didn’t check the 2nd source term since it should be the same
situation). I did this in kind of a hurry so maybe I am goofing up. I have
not studied the finite volume method, and am just going on a hunch.
It also seemed to me that there might be a minus sign missing in the
RobinCoeff. The solution I get is still pinned at initial condition
regardless of a minus sign or no minus sign there.
Thanks
On Fri, Dec 7, 2018 at 6:05 PM Guyer, Jonathan E. Dr. (Fed) via fipy <
fipy@nist.gov> wrote:
> Drew -
>
>
> Thanks to your notes, I found a couple of errors in our Robin discussion.
> As you noted, the units don't work for the conversion of the divergence of
> gamma grad phi to the sum over faces (midway on page 1 of your notes). The
> left hand side should be integrated over volume. That line is expressing
> the discretization of the divergence theorem.
>
> I also discovered that the .divergence operator in FiPy doesn't work
> reliably on scalars. There's no real reason it should, but I thought I had
> convinced myself that it did. As a result, everything in the Robin terms
> needs to be multiplied by the surface normal
>
> Those changes to the FiPy documentation are [here](
> https://github.com/usnistgov/fipy/pull/615) for now.
>
> Beyond addition that, there are some changes necessary to put your
> boundary condition into Robin form. The point of the Robin condition is
> that it ties the gradient of the field to the value of the field.
> So, g isn't h(T - T_inf) / (-k); it's just -h T_inf / (-k).
> Likewise, \vec{a} isn't zero. Rather, \hat{n}\cdot\vec{a} = h/(-k).
>
> Here are the changes I made to your script to reflect these changes:
>
> https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/compare/master...guyer:master
>
> I think I got your normalized temperature accounted for properly, and put
> the thermal diffusivity and convection coefficient in the proper places,
> but it's worth checking through.
>
> With these changes, the temperature field is now rotationally symmetric
> (it wasn't before, which is why I had to multiply by the normal before
> taking the divergence).
>
> Heat seems to be fluxing in from the outside, so I probably have the sign
> wrong.
>
> I don't have Octave, so I have no idea how this compares with your
> analytical solution.
>
> - Jon
>
> > On Dec 4, 2018, at 2:34 PM, Drew Davidson wrote:
> >
> > Hi Dr. Guyer,
> >
> > First I tried getting rid of the square brackets in
> ConvectionTestProb
Re: how to apply convection boundary condition in 2D?
Hi Dr. Guyer, First I tried getting rid of the square brackets in ConvectionTestProblem2D_01.py (commit ‘changed square brackets to parenthesis in convection BC, but get same…’), but results are the same (still wrong). Next: As you directed, I took a look at: https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions Firefox was showing me raw latex rather than human-readable equations at that web address, so I made a version I could read using the kile editor: https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsage.pdf That material is rather challenging for me. I took a stab at it, resulting in: commit: ‘made a first attempt at a formulation in view of suggestions by fipy …’ ConvectionTestProblem2D_01_2ndTry.py https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/fb47960548650c994eb0c6f990e0db297566e174/ConvectionTestProblem2D_01_2ndTry.py a few handwritten notes indicating what I was thinking: https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsageConvectionBC.pdf Here is a bit of the code: #ref: https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions #warning: avoid confusion between convectionCoeff in that fipy documentation, which refers to terms involving "a", and convectionCoeff here, which refers to a heat transfer convection coefficient at a boundary Gamma0=D_thermal Gamma = FaceVariable(mesh=mesh, value=Gamma0) mask=surfaceFaces Gamma.setValue(0., where=mask) dPf = FaceVariable(mesh=mesh, value=mesh._faceToCellDistanceRatio * mesh.cellDistanceVectors) Af = FaceVariable(mesh=mesh, value=mesh._faceAreas) #RobinCoeff = (mask * Gamma0 * Af / (dPf.dot(a) + b)).divergence #a is zero in our case b=1. RobinCoeff = (mask * Gamma0 * Af / b).divergence #a is zero in our case #eq = (TransientTerm() == DiffusionTerm(coeff=Gamma) + RobinCoeff * g - ImplicitSourceTerm(coeff=RobinCoeff * mesh.faceNormals.dot(a))) #a is zero in our case # g in this formulation is -convectionCoeff/k*var, where var=T-T_infinity eq = (TransientTerm() == DiffusionTerm(coeff=Gamma) + ImplicitSourceTerm(RobinCoeff * -convectionCoeff/k)) Now the solution variable remains stuck at the initial condition, as if the boundary condition is not being applied. I am sort of out of my depth at this point. I was guessing that an ImplicitSourceTerm was the right thing to do, since 'g' in the Robin condition depends on the solution variable. I did mess around a little in IPython seeing if any terms are coming out all zeros. Again, I put everything at https://github.com/cashTangoTangoCash/Convection2DFiPyExample01. Thanks On Mon, Dec 3, 2018 at 4:21 PM Guyer, Jonathan E. Dr. (Fed) via fipy < fipy@nist.gov> wrote: > Drew - > > Apologies for the delayed reply. > > There are a couple of issues here: > > `-convectionCoeff/k*(var.faceValue-T_infinity)` describes a FiPy > FaceVariable. > > `[-convectionCoeff/k*(var.faceValue-T_infinity)]` is a single element > Python list that contains a FiPy FaceVariable. > > Multiplying that list by other elements has rather unpredictable results. > > In short, get rid of the square brackets. We do our best to treat Python > lists like FiPy vector fields, but there's only so much we can do. A list > that holds a numpy array is not a vector field. > > > Beyond that, what you've described looks like a Robin boundary condition > to me. Our best recommendation for Robin conditions is covered at: > > > https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions > > Please don't hesitate to ask for more clarification if this doesn't get > you where you need. > > - Jon > > > > > On Nov 16, 2018, at 11:55 PM, Drew Davidson > wrote: > > > > Hello, > > > > I am stuck in how to correctly apply a simple convection boundary > condition in FiPy, in the context of simple transient heat conduction in 2D. > > > > Is this correct: > > > var.faceGrad.constrain([-convectionCoeff/k*(var.faceValue-T_infinity)]*mesh.faceNormals,where=surfaceFaces) > > > > I have a 2D mesh generated in gmsh. The convection boundary condition is > applied to a curved boundary. > > > > The code and results are found at: > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/tree/master > > > > The project appears to have a heap of files, but It’s a really just a > simple 2D problem with a comparison to an analytical solution. The basic > script is ConvectionTestProblem2D_01.py. Report20181116.pdf shows current > results, which are clearly wrong. > > > > Am I c
how to apply convection boundary condition in 2D?
Hello, I am stuck in how to correctly apply a simple convection boundary condition in FiPy, in the context of simple transient heat conduction in 2D. Is this correct: var.faceGrad.constrain([-convectionCoeff/k*(var.faceValue-T_infinity)]*mesh.faceNormals,where=surfaceFaces) I have a 2D mesh generated in gmsh. The convection boundary condition is applied to a curved boundary. The code and results are found at: https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/tree/master The project appears to have a heap of files, but It’s a really just a simple 2D problem with a comparison to an analytical solution. The basic script is ConvectionTestProblem2D_01.py. Report20181116.pdf shows current results, which are clearly wrong. Am I correctly applying the convection boundary condition in terms of the FiPy syntax/language for this 2D problem with a gmsh mesh and a curved boundary? Thanks ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
interpolation taking too long in an adaptive meshing effort (remeshing)
Hello, Short story: I am posting because before I tried this, I was unable to find on this mailing list a clear statement/warning that such an interpolation can run slow (eat up computer time). Is there a faster way to interpolate? Long story: I thought of an adaptive meshing scheme for examples.phase.anisotropy after reading: Ferreira, Alexandre Furtado, Leonardo de Olivé Ferreira, and Abner da Costa Assis. “Numerical Simulation of the Solidification of Pure Melt by a Phase-Field Model Using an Adaptive Computation Domain.” *Journal of the Brazilian Society of Mechanical Sciences and Engineering* 33, no. 2 (June 2011): 125–30. https://doi.org/10.1590/S1678-5878201100022. Their adaptive meshing scheme is quite simple and the paper is short. They do not interpolate, because they are writing their own code from the ground up. But in FiPy, it seems necessary to interpolate a solution from an old mesh to a new mesh after remeshing. I wrote a script and it runs for tiny problems, but seems useless because the interpolation in remeshing takes way too long (which defeats the purpose of ~Ferreira et al 2011 adaptive meshing), with interpolation time exploding as mesh cell count increases. I didn’t record the cell count exactly, but let’s say at ~5E5 cells, after 3.5 hours in a single interpolation, I gave up and hit ctrl-c. 2500 cells took less than 2 seconds, 2E4 cells 100 seconds, 6E4 cells 900 seconds. For 5E5 cells, my 4 core 2013 desktop PC shows multiple cores in use, but most cores at 50% load; using about 1/3 of the 8 GB RAM. Interpolation means transferring a FiPy solution from a first unstructured mesh of triangles to a second unstructured mesh of triangles (2D): newPhase=CellVariable(mesh=newMesh, value=oldPhase(newMesh.cellCenters,order=1), hasOld=True) new_dT=CellVariable(mesh=newMesh, value=old_dT(newMesh.cellCenters,order=1), hasOld=True) Mesh 2 always has more cells than mesh 1. Eventually mesh 2 is the mesh in examples.phase.anisotropy. Does anybody know a way to speed this interpolation up, within FiPy? I don’t believe I can use Grid2D meshes because my final domain of interest is not rectangular (examples.phase.anisotropy is just a beginner problem to play with). The meshes are unstructured sets of triangular cells and are created using gmsh via pygmsh/meshio. Mesh density varies with x,y (nonuniform meshes). Google seemed to lead me to various interpolation tools, but I do not see a way to take the data outside of FiPy to do the interpolation and then bring the data back into FiPy. I looked at the FiPy CellVariable source code but it is beyond me to mess around with it. It looked like it might run fast since it seems to operate on only the nearest neighbors, but maybe it is slowed down by having to do one cell at a time? Thanks ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Remeshing and code speedup
Hi Carsten, Have you looked at: https://www.ctcms.nist.gov/fipy/documentation/EFFICIENCY.html http://nbviewer.jupyter.org/github/wd15/fipy-efficiency/blob/master/notebooks/FiPy-IPython.ipynb https://www.mail-archive.com/fipy@nist.gov/msg03180.html thread I struggled to install weave and get --inline working. I didn't use it much yet, but my initial impression was the speedup was modest. Remeshing seems like a great possibility for making problems such as examples.phase.anisotropy run faster (use higher mesh density at the interface), but I am afraid that FiPy architecture (lazy evaluation etc) could interfere. I hope you share what you find with remeshing. Thanks On Tue, Jul 24, 2018 at 10:56 AM Carsten Langrock wrote: > Thanks for pointing out the performance order of the solvers. I’ll try to > get pySparse to work to compare it other solvers. It’s also good to know > that I shouldn’t give up on 2D just yet;-) > > Regards, > Carsten > > _ > *Dipl.-Phys. Carsten Langrock, Ph.D.* > > Senior Research Scientist > Edward L. Ginzton Laboratory, Rm. 202 > Stanford University > > 348 Via Pueblo Mall > 94305 Stanford, CA > > Tel. (650) 723-0464 > Fax (650) 723-2666 > > Ginzton Lab Shipping Address: > James and Anna Marie Spilker Engineering and Applied Sciences Building > 04-040 > 348 Via Pueblo Mall > 94305 Stanford, CA > _ > > On Jul 24, 2018, at 6:11 AM, Guyer, Jonathan E. Dr. (Fed) < > jonathan.gu...@nist.gov> wrote: > > FiPy still does not support remeshing. > > As Dario said, choice of solver can make a big difference. I've not used > PyAMG much, but PySparse is dramatically faster than SciPy. PyTrilinos is > slower than PySparse, but enables you to solve in parallel. > > I've also found that 2D problems solve much better than the 1D performance > would lead you to believe. There's just a lot of overhead in setting up the > problem and the Python communication with the lower-level libraries. > > > On Jul 23, 2018, at 6:44 PM, Carsten Langrock > wrote: > > Hi, > > Thanks for the help with getting FiPy running under Linux! I am trying to > re-create a 1D nonlinear diffusion problem for which we have C++ code that > uses the implicit Thomas algorithm based on > > J. Weickert, B. Romerny, M. Viergever, "Efficient and Reliable Schemes > for Nonlinear Diffusion Filtering”, IEEE transactions on Image Processing, > vol.7, N03, page 398, March 1998 > > I have been able to get results in FiPy that match this code very closely > which was a great start. Our C++ code uses a fixed number of spatial points > and a fixed time step, but re-meshes space to most efficiently use the size > of the array; it increases the spatial step size by 2 whenever the > concentration at a particular point reaches a set threshold. I tried > implementing this in FiPy as well, but haven’t had much luck so far. I saw > an old mailing-list entry from 2011 where a user was told that FiPy wasn’t > meant to do remeshing. Is that still the case? > > I’d imagine one would somehow need to update the Grid1D object with the > new ‘dx’, but since the CellVariable that holds the solution was > initialized with that mesh object, I am not sure that such a change would > propagate in a sensible fashion. I think I know how to map the value of the > CellVariable to account for the change in ‘dx’ by > > array_size = 2000 > phi.value = numpy.concatenate((phi.value[1:array_size/2:2], > numpy.zeros(1500))) > > for the case when the initial variable holds 2000 spatial points. Maybe > there’s a more elegant way, but I think this works in principle. > > Another question would be execution speed. Right now, even when not > plotting the intermediate solutions, it takes many seconds on a very > powerful computer to run a simple diffusion problem. I am probably doing > something really wrong. I wasn’t expecting the code to perform as well as > the C++ code, but I had hoped to come within an order of magnitude. Are > there ways to optimize the performance? Maybe select a particularly clever > solver? If someone could point me into the right direction that’d be great. > In the end, I would like to expand the code into 2D, but given the poor 1D > performance, I don’t think that this would be feasible at this point. > > Thanks, > Carsten > > _ > Dipl.-Phys. Carsten Langrock, Ph.D. > > Senior Research Scientist > Edward L. Ginzton Laboratory, Rm. 202 > Stanford University > > 348 Via Pueblo Mall > 94305 Stanford, CA > > Tel. (650) 723-0464 > Fax (650) 723-2666 > > Ginzton Lab Shipping Address: > James and Anna Marie Spilker Engineering and Applied Sciences Building > 04-040 > 348 Via Pueblo Mall > 94305 Stanford, CA > _ > ___ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > > ___
Re: thermal contact between two regions
[0] ==
L_Substrate))
v = fp.CellVariable(mesh=mesh)
v.constrain(0, where=mesh.facesLeft)
v.faceGrad.constrain(HeatFlux / k_Sample, where=mesh.facesRight)
#follow examples at
https://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html#module-examples.diffusion.mesh1D
eq1=fp.TransientTerm()==fp.DiffusionTerm(coeff=diffCoeff) #implicit
dtExplicitLimit_Substrate = cellSize**2 / (2 * D_Substrate)
dtExplicitLimit_Sample = cellSize**2 / (2 * D_Sample)
dtExplicitLimit_InterfaceCell=cellSize**2/(2*Keff.value.mean()/rho_times_cp_eff)
#play with Keff in ipython; .mean seems close enough to what you want
dt=10*min([dtExplicitLimit_Substrate,dtExplicitLimit_Sample,dtExplicitLimit_InterfaceCell])
#SETTING
#viewer=fp.Viewer(vars=(v,),datamin=0,datamax=10.)
viewer=fp.Viewer(vars=(v,))
TJumpAtInterfaceExpected=HeatFlux*R_TIM_unit_area
TJumpInBarNoInterfaceExpected=HeatFlux*(L_Substrate/k_Substrate+L_Sample/k_Sample)
TJumpInBarExpected=TJumpAtInterfaceExpected+TJumpInBarNoInterfaceExpected
blurb='''
at steady state:
expected T jump at interface is %.2E
expected T jump in bar with no interface is %.2E
expected total T jump in bar is %.2E
hit return to continue
''' %
(TJumpAtInterfaceExpected,TJumpInBarNoInterfaceExpected,TJumpInBarExpected)
#run with --inline if you have weave working
steps = 15000 #SETTING
plotEvery=100 #SETTING
saveDataEvery=200 #SETTING
#TODO need code to stop when at thermal steady state
times=[]
rod_dTs=[]
outFilenameBase='20180717UniformMesh04' #SETTING
outFilename=outFilenameBase+'.csv'
assert (not os.path.exists(outFilename)), '%s already exists on disk;
this script requires that you manually update outFilenameBase'
with open(outFilename,'w') as csvfile:
csvWriter = csv.writer(csvfile, delimiter=' ',quotechar='|',
quoting=csv.QUOTE_MINIMAL)
for step in range(steps):
eq1.solve(var=v,dt=dt)
if step%plotEvery==0:
print 'step is %d of %d\n' % (step,steps)
viewer.plot()
if step%saveDataEvery==0:
#saving the temperature drop across the entire rod vs time:
times.append(step*dt)
rod_dTs.append(v.value[-1]-v.value[0])
csvWriter.writerow([step*dt,v.value[-1]-v.value[0]]) #TODO
am looking at values at cell centers; really want face values
raw_input(blurb)
fig, ax = plt.subplots()
ax.plot(times,rod_dTs)
ax.set(xlabel='time (s)', ylabel='temperature drop across entire rod,
K')
ax.grid()
fig.savefig(outFilenameBase+'.png')
plt.show()
raw_input('hit enter to continue')
#make a copy of this script with a new name so you have a good record
of what produced the result
#copyfile(src,dst)
copyfile('ModifyWheelerInViewOfGuyer20130517TransientVersionImplicit03.py',outFilenameBase+'.py')
#don't forget to run with --inline if you have weave installed and
working
On Mon, Jul 23, 2018 at 11:30 AM Daniel Wheeler
wrote:
> On Tue, Jul 17, 2018 at 4:43 PM, Drew Davidson
> wrote:
> >
> > I still need FiPy code for:
> >
> >
> > dAP1/(dAP1+dAP2)
> >
> >
> > and
> >
> >
> > dAP2/(dAP1+dAP2)
> >
> >
> > where dAP1 and dAP2 were distances from cell center to cell face for
> cells
> > on either side of the interface. If these FiPy expressions are
> unavailable,
> > I would think assuming .5 is going to be OK...
>
> I think it's outlined in the link that you gave.
>
> Kcontact = hc * mesh._cellDistances
> Kavg = Kcell.harmonicFaceValue
> K.setValue(Kcontact * Kavg / (Kavg + Kcontact * (1 - dx /
> mesh._cellDistances)), where=mesh.physicalFaces['thermal contact'])
>
> m._cellDistances is the distance between each cell and `dx` is the
> volume of the cell for a 1D problem.
>
>
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
how to compute the area of the solid material in examples.phase.anisotropy?
Hello, Would anyone please give the FiPy code which computes the area of the solidified material in examples.phase.anisotropy? Searching for --area-- in the Mailing List Archives returns a rather broad result set. Thanks ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: thermal contact between two regions
Hello, I continued to stare at https://www.mail-archive.com/fipy@nist.gov/msg02641.html. I realized Keff must be effective thermal conductivity as described in: Salazar, Agust n. “On Thermal Diffusivity.” *European Journal of Physics* 24, no. 4 (July 1, 2003): 351–58. https://doi.org/10.1088/0143-0807/24/4/353 . For years I sat right next to grad students whose entire degree was in computing Keff of composite materials via ANSYS/COMSOL, but I still looked at https://www.mail-archive.com/fipy@nist.gov/msg02641.html with incomprehension. I realized all that was being done in https://www.mail-archive.com/fipy@nist.gov/msg02641.html must be to compute Keff of the finite volume which contains the thermal contact resistance, and then to assign thermal conductivity of the finite volume face at the interface to Keff. The Salazar paper also gives equations for effective thermal diffusivity, which is what I was asking for in the current thread. Now I think I mostly have what I need. I still need FiPy code for: dAP1/(dAP1+dAP2) and dAP2/(dAP1+dAP2) where dAP1 and dAP2 were distances from cell center to cell face for cells on either side of the interface. If these FiPy expressions are unavailable, I would think assuming .5 is going to be OK... Thanks On Thu, Jul 12, 2018 at 2:28 PM Daniel Wheeler wrote: > I think you can do this by aligning the mesh with the contact region > and specifying the diffusion coefficient at the contact faces. > > On Thu, Jul 12, 2018 at 12:28 PM, Drew Davidson > wrote: > > Hello, > > > > > > I’m sorry for not being more specific. I was trying to incorporate > > https://www.mail-archive.com/fipy@nist.gov/msg02626.html by reference. > > -- > Daniel Wheeler > > ___ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Can Gmsh2D read a 1D gmsh mesh?
Hello,
It seems possible to generate a parameterized 1D mesh using gmsh within a
python script, and then generate the corresponding 1D FiPy mesh in that
same python script:
1) Use pygmsh in a python script to construct a parameterized 1D (gmsh)
mesh.
2) From that 1D mesh, construct a list of dx's
3) Follow https://www.mail-archive.com/fipy@nist.gov/msg01650.html and
create a 1D FiPy mesh from a list of dx's.
It was pure luck to have found the mailing list message in (3).
I did not see how to get the Physical Groups data from the pygmsh mesh into
the FiPy mesh.
Thanks
On Thu, Jul 12, 2018 at 3:26 PM Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:
> No, FiPy doesn't presently have any ability to read 1D meshes.
>
> I do have some experimental code that supports this, but it needs a lot of
> work. I have no prognosis for when this will be made public.
>
> - Jon
>
> > On Jul 12, 2018, at 1:21 PM, Drew Davidson
> wrote:
> >
> > Hello,
> >
> > Can Gmsh2D read a 1D gmsh mesh?
> >
> > I have tried a similar script to the FiPy diffusion example circle.py:
> >
> > SCRIPT BEGINS
> > import fipy as fp
> > from IPython import embed
> > import sys
> >
> > #I created this gmsh .geo file in gmsh GUI, then copied and pasted here;
> I cannot see a syntax error if it exists???
> >
> > geo="""
> > L_Sample=%(L_Sample)g;
> > L_Substrate=%(L_Substrate)g;
> > cellSize=%(cellSize)g;
> > Point(1) = {0, 0, 0, cellSize};
> > Point(2) = {L_Substrate, 0, 0, cellSize};
> > Point(3) = {L_Substrate+L_Sample, 0, 0, cellSize};
> > Line(1) = {1, 2};
> > Line(2) = {2, 3};
> > Physical Point("SubstrateEnd") = {1};
> > Physical Point("SampleEnd") = {3};
> > Physical Point("ThermalContact") = {2};
> > Physical Line("Substrate") = {1};
> > Physical Line("Sample") = {2};
> > """ % locals()
> >
> > #I am not having success with:
> > mesh = fp.Gmsh2D(geo,coordDimensions=1) #TODO throwing error about there
> being no cells; is it a problem that mesh is 1D and Gmsh2D is called 2D?
> > mesh = fp.Gmsh2D(geo) #TODO throwing error about there being no cells
> >
> > #then I went into gmsh GUI and saved to .msh, but still no luck:
> > mesh = fp.Gmsh2D('OneDimensionalSampleAndSubstrate.msh') #TODO throwing
> error about there being no cells
> > mesh =
> fp.Gmsh2D('OneDimensionalSampleAndSubstrate.msh',coordDimensions=1) #TODO
> throwing error about there being no cells
> > SCRIPT ENDS
> >
> >
> > If coordDimensions=1 is unsupported or meaningless, it would be great if
> FiPy printed a message saying so.
> >
> > I believe my FiPy installation is OK and was downloaded from github
> 6/30/2018 (fipy-develop). I can run circle.py, although I have never gotten
> it to plot the mesh. I can run anisotropy.py from diffusivity examples.
> Nearly all fipy tests pass (setup.py). Can you please check if Gmsh2D is
> capable of reading a 1D gmsh mesh? Or is there another FiPy class which can
> do it, and how?
> >
> > gmsh –version returns 3.0.6. Installation is system-wide via download
> method in Ubuntu 16.04.
> >
> > If I cannot read in a parameterized 1D gmsh mesh using FiPy
> functionality, I am at a loss as to how to have a single Python script that
> allows me to set cellSize, L_Sample, and L_Substrate, produce a mesh, and
> efficiently conduct a parameter study with FiPy. Gmsh is desirable because
> I want to go to a more complicated nonuniform/graded mesh, and would rather
> not write code to make a mesh myself.
> >
> > Thanks.
> >
> >
> > ___
> > fipy mailing list
> > fipy@nist.gov
> > http://www.ctcms.nist.gov/fipy
> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Can Gmsh2D read a 1D gmsh mesh?
Hello,
Can Gmsh2D read a 1D gmsh mesh?
I have tried a similar script to the FiPy diffusion example circle.py:
SCRIPT BEGINS
import fipy as fp
from IPython import embed
import sys
#I created this gmsh .geo file in gmsh GUI, then copied and pasted here; I
cannot see a syntax error if it exists???
geo="""
L_Sample=%(L_Sample)g;
L_Substrate=%(L_Substrate)g;
cellSize=%(cellSize)g;
Point(1) = {0, 0, 0, cellSize};
Point(2) = {L_Substrate, 0, 0, cellSize};
Point(3) = {L_Substrate+L_Sample, 0, 0, cellSize};
Line(1) = {1, 2};
Line(2) = {2, 3};
Physical Point("SubstrateEnd") = {1};
Physical Point("SampleEnd") = {3};
Physical Point("ThermalContact") = {2};
Physical Line("Substrate") = {1};
Physical Line("Sample") = {2};
""" % locals()
#I am not having success with:
mesh = fp.Gmsh2D(geo,coordDimensions=1) #TODO throwing error about there
being no cells; is it a problem that mesh is 1D and Gmsh2D is called 2D?
mesh = fp.Gmsh2D(geo) #TODO throwing error about there being no cells
#then I went into gmsh GUI and saved to .msh, but still no luck:
mesh = fp.Gmsh2D('OneDimensionalSampleAndSubstrate.msh') #TODO throwing
error about there being no cells
mesh = fp.Gmsh2D('OneDimensionalSampleAndSubstrate.msh',coordDimensions=1)
#TODO throwing error about there being no cells
SCRIPT ENDS
If coordDimensions=1 is unsupported or meaningless, it would be great if
FiPy printed a message saying so.
I believe my FiPy installation is OK and was downloaded from github
6/30/2018 (fipy-develop). I can run circle.py, although I have never gotten
it to plot the mesh. I can run anisotropy.py from diffusivity examples.
Nearly all fipy tests pass (setup.py). Can you please check if Gmsh2D is
capable of reading a 1D gmsh mesh? Or is there another FiPy class which can
do it, and how?
gmsh –version returns 3.0.6. Installation is system-wide via download
method in Ubuntu 16.04.
If I cannot read in a parameterized 1D gmsh mesh using FiPy functionality,
I am at a loss as to how to have a single Python script that allows me to
set cellSize, L_Sample, and L_Substrate, produce a mesh, and efficiently
conduct a parameter study with FiPy. Gmsh is desirable because I want to go
to a more complicated nonuniform/graded mesh, and would rather not write
code to make a mesh myself.
Thanks.
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: thermal contact between two regions
Hello, I’m sorry for not being more specific. I was trying to incorporate https://www.mail-archive.com/fipy@nist.gov/msg02626.html by reference. Sketch of problem: (1) Region 1 --* thermal contact resistance * -- Region 2 (2) Temperature at (1) is zero. Heat flux at (2) is constant. At time zero, temperature is everywhere zero. ** What is a thermal contact resistance? https://en.wikipedia.org/wiki/Thermal_contact_conductance A dissertation available free online: Thermal contact resistance Author: Mikic, B. B.; Rohsenow, Warren M. Citable URI: http://hdl.handle.net/1721.1/61449 It’s equation (1.1) on page 12 of this dissertation. Thermal contact resistance is directly analogous to an electrical resistor in a simple circuit: (temperature drop across thermal contact resistance) = (heat flux) * (its unit-area thermal resistance) Units: Kelvin = W/m^2 times Kelvin/(W/m^2) Mathematically, the thermal contact resistance has no thickness and no heat capacity (it does not store energy). ** Please look at https://www.mail-archive.com/fipy@nist.gov/msg02641.html. Dr. Guyer and Dr. Wheeler considered the case of steady heat flow (steady state problem with no transient term), and their thermal contact has a finite thickness ‘dx’. They derived an expression for effective thermal conductivity at the thermal contact location, and provided FiPy code. How do things change in a transient problem? What is the effective thermal diffusivity (thermal diffusivity = thermal conductivity/(mass density * specific heat)) at the interface location, given that Region 1 and Region 2 have different thermal conductivities, mass densities, and specific heats? How is it written in FiPy? This is the point at which my knowledge of the finite volume method and FiPy is too weak. A solution where the thermal contact has a finite thickness ‘dx’ is also of interest, since it seems like a thermal contact resistance can be approximated by setting ‘dx’ very small. Trying to get on the same page: thermal diffusivity m^2/s thermal conductivity W/m/K mass density kg/m^3 specific heat J/kg/K Thanks On Thu, Jul 12, 2018 at 9:49 AM Daniel Wheeler wrote: > On Wed, Jul 11, 2018 at 7:00 PM, Drew Davidson > wrote: > > > > 1. Transient heat conduction rather than steady state. Same boundary > > conditions but initial condition can be zero temperature eveywhere. > > Add a TransientTerm to the equation and step through the solution with > time steps and sweeps if necessary. > > > 2. The two regions have differing properties (thermal conductivity, mass > > density, specific heat) > > The coefficients can have spatially varying values in FiPy. Define the > coefficients as face or cell variables. > > > 3. The thermal contact is a true thermal contact resistance (a pure > > resistance); it has zero thickness, and it has zero heat capacity. > > I need to see it defined mathematically to understand how to implement > it in FiPy, but I'm confident it's possible to define it. > > -- > Daniel Wheeler > ___ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
thermal contact between two regions
Hello, I want to do a heat conduction problem with FiPy that is similar to that in https://www.mail-archive.com/fipy@nist.gov/msg02626.html “how to implement a thermal contact between two regions?” What is different from that thread: 1. Transient heat conduction rather than steady state. Same boundary conditions but initial condition can be zero temperature eveywhere. 2. The two regions have differing properties (thermal conductivity, mass density, specific heat) 3. The thermal contact is a true thermal contact resistance (a pure resistance); it has zero thickness, and it has zero heat capacity. I tried to guess at code from looking at the aforementioned mailing list thread, but my understanding of FiPy is too weak; it would be great to see an authoritative answer. Thanks ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to compute interface velocity in examples.phase.anisotropy?
Hello,
I was just looking at:
Boettinger, W. J., J. A. Warren, C. Beckermann, and A. Karma. “Phase-Field
Simulation of Solidification.” Annual Review of Materials Research 32, no.
1 (August 1, 2002): 163–94.
https://doi.org/10.1146/annurev.matsci.32.101901.155803.
and maybe Equation 27 is the way to compute interface velocity in
examples.phase.anisotropy:
v = -\frac{1}{\| \nabla \phi \|}\frac{\partial \phi}{\partial t}
Then it remains to express the right hand side of that equation in the FiPy
language. It does look like this only gives the velocity magnitude rather
than the velocity vector.
On Tue, Apr 3, 2018 at 12:27 PM, Drew Davidson
wrote:
> Hello,
>
> In FiPy, how would one calculate the interface velocity at a given time
> step in examples.phase.anisotropy?
>
> At this time, I was mostly interested in the max and min of the magnitude
> of the interface velocity.
>
> Thanks
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
how to compute interface velocity in examples.phase.anisotropy?
Hello, In FiPy, how would one calculate the interface velocity at a given time step in examples.phase.anisotropy? At this time, I was mostly interested in the max and min of the magnitude of the interface velocity. Thanks ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
should the links on this page work?
Hi, I came across this page: https://github.com/usnistgov/fipy/wiki/FiPyExamples The links seem to be described as an additional set of FiPy instructional materials, perhaps authored by non-NIST university personnel. When I click these links, I get a Resource Not Found message in a tab with a NIST logo. I didn't click all links, but I clicked quite a few and always got Resource Not Found. What is the status of this material? Is the general public able to access it at this time, and how? Thanks. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Fwd: axisymmetric solidification problem in FiPy: trying to start with a basic 1D Stefan problem (solidification)
-- Forwarded message -- From: Drew Davidson Date: Mon, Feb 13, 2017 at 12:33 PM Subject: axisymmetric solidification problem in FiPy: trying to start with a basic 1D Stefan problem (solidification) To: fipy@nist.gov Hi, Sorry for the two emails; it's unclear whether the address is fipy@nist.gov or mailto:fipy@nist.gov. 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="t
axisymmetric solidification problem in FiPy: trying to start with a basic 1D Stefan problem (solidification)
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 vel
