ANN: FiPy 3.2
We are pleased to announce the release of FiPy 3.2. http://www.ctcms.nist.gov/fipy This is predominantly a DevOps release. The focus has been on making FiPy easier to install with conda. It's also possible to install a minimal set of prerequisites with pip. Further, FiPy is automatically tested on all major platforms using cloud-based Continuous Integration (linux with CircleCI, macOS with TravisCI, and Windows with AppVeyor). This release addresses 38 issues. FiPy is an object oriented, partial differential equation (PDE) solver, written in Python, based on a standard finite volume (FV) approach. The framework has been developed in the Metallurgy Division and Center for Theoretical and Computational Materials Science (CTCMS), in the Material Measurement Laboratory (MML) at the National Institute of Standards and Technology (NIST). The solution of coupled sets of PDEs is ubiquitous to the numerical simulation of science problems. Numerous PDE solvers exist, using a variety of languages and numerical approaches. Many are proprietary, expensive and difficult to customize. As a result, scientists spend considerable resources repeatedly developing limited tools for specific problems. Our approach, combining the FV method and Python, provides a tool that is extensible, powerful and freely available. A significant advantage to Python is the existing suite of tools for array calculations, sparse matrices and data rendering. The FiPy framework includes terms for transient diffusion, convection and standard sources, enabling the solution of arbitrary combinations of coupled elliptic, hyperbolic and parabolic PDEs. Currently implemented models include phase field treatments of polycrystalline, dendritic, and electrochemical phase transformations as well as a level set treatment of the electrodeposition process. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: 1D diffusion problem
Hi Jon,
Thank you so much for taking the time to explain everything!
Clara
> On 22 Apr 2019, at 13:31, Guyer, Jonathan E. Dr. (Fed) via fipy
> wrote:
>
> Also, the natural boundary condition for FiPy is no-flux, so there's no need
> to do
>
> T_var.faceGrad.constrain(0, where=mesh.facesLeft)
>
>> On Apr 22, 2019, at 1:28 PM, Guyer, Jonathan E. Dr. (Fed) via fipy
>> wrote:
>>
>> Clara -
>>
>> - The cause of everything dropping to zero is that T_var.old is initialized
>> to zero. This isn't a particularly useful initial value, but that's what it
>> is. Usually, T_var.updateOld() is the first thing called in the time step
>> loop, so this initial condition doesn't matter, but because of the way
>> you've structured your adaptive time steps, the first sweep is done with an
>> erroneous initial condition.
>>
>> Call T_var.updateOld() one time before you start taking time steps.
>>
>> - Your sweep tolerance is too severe. I adjusted to 1.0e-6 and the adaptive
>> stepper works much better.
>>
>> - The way you define and update dt, you'll get recursion errors eventually.
>> Instead, define
>>
>> dt = 0.9*dr**2/(2.0*Diff[0].value)
>>
>> If you're curious, this is because:
>>
>> * Diff is a FaceVariable
>> * so, Diff[0] is a UnaryOperator Variable
>> * so, dt = 0.9*dr**2/(2.0*Diff[0]) is a BinaryOperator Variable
>> * so, dt *= 1.1/dt *= 0.8 is a BinaryOperator Variable that gets nested
>> deeper
>>and deeper every time step
>>
>> Diff[0].value short circuits that
>>
>> - Your problem is completely linear, so sweeps aren't necessary at this point
>>
>> - This equation is fully implicit, so there's no need to start with such a
>> tiny initial time step
>>
>> - Jon
>>
>>> On Apr 19, 2019, at 2:42 PM, Clara Maurel wrote:
>>>
>>> Hello everyone,
>>>
>>> I apologize in advance for this very basic question…
>>> I am trying to solve a 1D diffusion problem where I have a cooling planet
>>> (core+mantle) with an initial temperature profile and a diffusion
>>> coefficient that takes different values in the core and in the mantle. I
>>> want a fixed boundary condition at the surface (T = 200 K) and a zero
>>> gradient at the center of the body, but the temperature at the center (left
>>> of the mesh) is not fixed and should cool.
>>> When I run the script, the whole profile immediately goes down to zero
>>> (except the end right of the mesh) and diffuses from there. I must be doing
>>> something wrong with the left boundary condition but I can’t figure out
>>> what.
>>>
>>> Could someone help me to find out what I missed? Any hint would be greatly
>>> appreciated!
>>>
>>> Thank you in advance for your help!
>>> Clara
>>>
>>> Here is my simple script:
>>>
>>> from datetime import datetime
>>> startTime = datetime.now()
>>> import numpy as np
>>> import scipy
>>> from fipy import *
>>> from fipy.solvers.pysparse import LinearLUSolver as Solver
>>> from sympy import *
>>> from scipy import interpolate
>>>
>>> def Get_closest_id(L,value):
>>> return list(L).index(min(L, key=lambda x:abs(x-value)))
>>>
>>> ##
>>> # DIFFUSION EQUATION: INITIALIZATION #
>>> ##
>>>
>>> ## Dimensions (km)
>>> R_body = 200.0
>>> R_core = 50.0
>>>
>>> ## Temperatures (K)
>>> T_core = 1200.0
>>> T_surf = 200.0
>>>
>>> ## Mesh
>>> nr = 1000
>>> dr = R_body / nr
>>> mesh = Grid1D(nx=nr, dx=dr)
>>>
>>> ## Variable
>>> T_var = CellVariable(name="T", mesh=mesh, hasOld=True)
>>>
>>> ## Initial temperature profile
>>> slope = (T_core-T_surf)/(R_core-R_body)
>>> intercept = T_core-slope*R_core
>>> T_var[:] = np.array([T_core]*int(R_core/dr) +
>>> list(slope*np.arange(R_core+dr, R_body+dr, dr)+intercept))
>>>
>>> ## Initial conditions (fixed surface value, null gradient at the center of
>>> the core)
>>> T_var.constrain(T_surf, mesh.facesRight)
>>> T_var.faceGrad.constrain(0, where=mesh.facesLeft)
>>>
>>> ## Diffusion coefficient (different in core and mantle, in km2 My-1)
>>> Diff = FaceVariable(mesh=mesh)
>>> X = mesh.faceCenters[0]
>>> Diff.setValue(150.0, where=(R_core >= X))
>>> Diff.setValue(15.0, where=(R_core < X) & (R_body >= X))
>>>
>>> ## Equation
>>> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff)
>>>
>>> if __name__ == "__main__":
>>> viewer = Viewer(vars=(T_var),ymin=0,ymax=1400)
>>> viewer.plot()
>>>
>>> if __name__ == '__main__':
>>> raw_input("Press to proceed...")
>>>
>>> ## Time and integration parameters (time in My; 3.154e13 = 1 My in s)
>>> time = 0.0
>>> duration = 1000.0
>>> dt = 0.9*dr**2/(2.0*Diff[0])
>>>
>>> total_sweeps = 4
>>> tolerance = 1.0e-10
>>>
>>> solver = Solver()
>>>
>>> #---#
>>> # DIFFUSION EQUATION: MAIN LOOP #
>>> #---#
>>>
>>> while time < duration:
>>>
>>> res0 = eq.sweep(T_var, dt=dt, solver=solver)
>>>
>>> for sweeps in range(total_sweeps):
>>> res = eq.sweep(T_var, dt=dt, solver=solver)
Re: 1D diffusion problem
Also, the natural boundary condition for FiPy is no-flux, so there's no need to
do
T_var.faceGrad.constrain(0, where=mesh.facesLeft)
> On Apr 22, 2019, at 1:28 PM, Guyer, Jonathan E. Dr. (Fed) via fipy
> wrote:
>
> Clara -
>
> - The cause of everything dropping to zero is that T_var.old is initialized
> to zero. This isn't a particularly useful initial value, but that's what it
> is. Usually, T_var.updateOld() is the first thing called in the time step
> loop, so this initial condition doesn't matter, but because of the way you've
> structured your adaptive time steps, the first sweep is done with an
> erroneous initial condition.
>
> Call T_var.updateOld() one time before you start taking time steps.
>
> - Your sweep tolerance is too severe. I adjusted to 1.0e-6 and the adaptive
> stepper works much better.
>
> - The way you define and update dt, you'll get recursion errors eventually.
> Instead, define
>
> dt = 0.9*dr**2/(2.0*Diff[0].value)
>
> If you're curious, this is because:
>
> * Diff is a FaceVariable
> * so, Diff[0] is a UnaryOperator Variable
> * so, dt = 0.9*dr**2/(2.0*Diff[0]) is a BinaryOperator Variable
> * so, dt *= 1.1/dt *= 0.8 is a BinaryOperator Variable that gets nested
> deeper
> and deeper every time step
>
> Diff[0].value short circuits that
>
> - Your problem is completely linear, so sweeps aren't necessary at this point
>
> - This equation is fully implicit, so there's no need to start with such a
> tiny initial time step
>
> - Jon
>
>> On Apr 19, 2019, at 2:42 PM, Clara Maurel wrote:
>>
>> Hello everyone,
>>
>> I apologize in advance for this very basic question…
>> I am trying to solve a 1D diffusion problem where I have a cooling planet
>> (core+mantle) with an initial temperature profile and a diffusion
>> coefficient that takes different values in the core and in the mantle. I
>> want a fixed boundary condition at the surface (T = 200 K) and a zero
>> gradient at the center of the body, but the temperature at the center (left
>> of the mesh) is not fixed and should cool.
>> When I run the script, the whole profile immediately goes down to zero
>> (except the end right of the mesh) and diffuses from there. I must be doing
>> something wrong with the left boundary condition but I can’t figure out what.
>>
>> Could someone help me to find out what I missed? Any hint would be greatly
>> appreciated!
>>
>> Thank you in advance for your help!
>> Clara
>>
>> Here is my simple script:
>>
>> from datetime import datetime
>> startTime = datetime.now()
>> import numpy as np
>> import scipy
>> from fipy import *
>> from fipy.solvers.pysparse import LinearLUSolver as Solver
>> from sympy import *
>> from scipy import interpolate
>>
>> def Get_closest_id(L,value):
>>return list(L).index(min(L, key=lambda x:abs(x-value)))
>>
>> ##
>> # DIFFUSION EQUATION: INITIALIZATION #
>> ##
>>
>> ## Dimensions (km)
>> R_body = 200.0
>> R_core = 50.0
>>
>> ## Temperatures (K)
>> T_core = 1200.0
>> T_surf = 200.0
>>
>> ## Mesh
>> nr = 1000
>> dr = R_body / nr
>> mesh = Grid1D(nx=nr, dx=dr)
>>
>> ## Variable
>> T_var = CellVariable(name="T", mesh=mesh, hasOld=True)
>>
>> ## Initial temperature profile
>> slope = (T_core-T_surf)/(R_core-R_body)
>> intercept = T_core-slope*R_core
>> T_var[:] = np.array([T_core]*int(R_core/dr) +
>> list(slope*np.arange(R_core+dr, R_body+dr, dr)+intercept))
>>
>> ## Initial conditions (fixed surface value, null gradient at the center of
>> the core)
>> T_var.constrain(T_surf, mesh.facesRight)
>> T_var.faceGrad.constrain(0, where=mesh.facesLeft)
>>
>> ## Diffusion coefficient (different in core and mantle, in km2 My-1)
>> Diff = FaceVariable(mesh=mesh)
>> X = mesh.faceCenters[0]
>> Diff.setValue(150.0, where=(R_core >= X))
>> Diff.setValue(15.0, where=(R_core < X) & (R_body >= X))
>>
>> ## Equation
>> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff)
>>
>> if __name__ == "__main__":
>>viewer = Viewer(vars=(T_var),ymin=0,ymax=1400)
>>viewer.plot()
>>
>> if __name__ == '__main__':
>>raw_input("Press to proceed...")
>>
>> ## Time and integration parameters (time in My; 3.154e13 = 1 My in s)
>> time = 0.0
>> duration = 1000.0
>> dt = 0.9*dr**2/(2.0*Diff[0])
>>
>> total_sweeps = 4
>> tolerance = 1.0e-10
>>
>> solver = Solver()
>>
>> #---#
>> # DIFFUSION EQUATION: MAIN LOOP #
>> #---#
>>
>> while time < duration:
>>
>>res0 = eq.sweep(T_var, dt=dt, solver=solver)
>>
>>for sweeps in range(total_sweeps):
>>res = eq.sweep(T_var, dt=dt, solver=solver)
>>
>>if res < res0 * tolerance:
>>time += dt
>>dt *= 1.1
>>T_var.updateOld()
>>print dt, res
>>if __name__ == '__main__':
>>viewer.plot()
>>
>>else:
>>dt *= 0.8
>>T_var[:] = T_var.old
>>print "warning " + str(res)
Re: 1D diffusion problem
Clara -
- The cause of everything dropping to zero is that T_var.old is initialized to
zero. This isn't a particularly useful initial value, but that's what it is.
Usually, T_var.updateOld() is the first thing called in the time step loop, so
this initial condition doesn't matter, but because of the way you've structured
your adaptive time steps, the first sweep is done with an erroneous initial
condition.
Call T_var.updateOld() one time before you start taking time steps.
- Your sweep tolerance is too severe. I adjusted to 1.0e-6 and the adaptive
stepper works much better.
- The way you define and update dt, you'll get recursion errors eventually.
Instead, define
dt = 0.9*dr**2/(2.0*Diff[0].value)
If you're curious, this is because:
* Diff is a FaceVariable
* so, Diff[0] is a UnaryOperator Variable
* so, dt = 0.9*dr**2/(2.0*Diff[0]) is a BinaryOperator Variable
* so, dt *= 1.1/dt *= 0.8 is a BinaryOperator Variable that gets nested
deeper
and deeper every time step
Diff[0].value short circuits that
- Your problem is completely linear, so sweeps aren't necessary at this point
- This equation is fully implicit, so there's no need to start with such a tiny
initial time step
- Jon
> On Apr 19, 2019, at 2:42 PM, Clara Maurel wrote:
>
> Hello everyone,
>
> I apologize in advance for this very basic question…
> I am trying to solve a 1D diffusion problem where I have a cooling planet
> (core+mantle) with an initial temperature profile and a diffusion coefficient
> that takes different values in the core and in the mantle. I want a fixed
> boundary condition at the surface (T = 200 K) and a zero gradient at the
> center of the body, but the temperature at the center (left of the mesh) is
> not fixed and should cool.
> When I run the script, the whole profile immediately goes down to zero
> (except the end right of the mesh) and diffuses from there. I must be doing
> something wrong with the left boundary condition but I can’t figure out what.
>
> Could someone help me to find out what I missed? Any hint would be greatly
> appreciated!
>
> Thank you in advance for your help!
> Clara
>
> Here is my simple script:
>
> from datetime import datetime
> startTime = datetime.now()
> import numpy as np
> import scipy
> from fipy import *
> from fipy.solvers.pysparse import LinearLUSolver as Solver
> from sympy import *
> from scipy import interpolate
>
> def Get_closest_id(L,value):
> return list(L).index(min(L, key=lambda x:abs(x-value)))
>
> ##
> # DIFFUSION EQUATION: INITIALIZATION #
> ##
>
> ## Dimensions (km)
> R_body = 200.0
> R_core = 50.0
>
> ## Temperatures (K)
> T_core = 1200.0
> T_surf = 200.0
>
> ## Mesh
> nr = 1000
> dr = R_body / nr
> mesh = Grid1D(nx=nr, dx=dr)
>
> ## Variable
> T_var = CellVariable(name="T", mesh=mesh, hasOld=True)
>
> ## Initial temperature profile
> slope = (T_core-T_surf)/(R_core-R_body)
> intercept = T_core-slope*R_core
> T_var[:] = np.array([T_core]*int(R_core/dr) + list(slope*np.arange(R_core+dr,
> R_body+dr, dr)+intercept))
>
> ## Initial conditions (fixed surface value, null gradient at the center of
> the core)
> T_var.constrain(T_surf, mesh.facesRight)
> T_var.faceGrad.constrain(0, where=mesh.facesLeft)
>
> ## Diffusion coefficient (different in core and mantle, in km2 My-1)
> Diff = FaceVariable(mesh=mesh)
> X = mesh.faceCenters[0]
> Diff.setValue(150.0, where=(R_core >= X))
> Diff.setValue(15.0, where=(R_core < X) & (R_body >= X))
>
> ## Equation
> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff)
>
> if __name__ == "__main__":
> viewer = Viewer(vars=(T_var),ymin=0,ymax=1400)
> viewer.plot()
>
> if __name__ == '__main__':
> raw_input("Press to proceed...")
>
> ## Time and integration parameters (time in My; 3.154e13 = 1 My in s)
> time = 0.0
> duration = 1000.0
> dt = 0.9*dr**2/(2.0*Diff[0])
>
> total_sweeps = 4
> tolerance = 1.0e-10
>
> solver = Solver()
>
> #---#
> # DIFFUSION EQUATION: MAIN LOOP #
> #---#
>
> while time < duration:
>
> res0 = eq.sweep(T_var, dt=dt, solver=solver)
>
> for sweeps in range(total_sweeps):
> res = eq.sweep(T_var, dt=dt, solver=solver)
>
> if res < res0 * tolerance:
> time += dt
> dt *= 1.1
> T_var.updateOld()
> print dt, res
> if __name__ == '__main__':
> viewer.plot()
>
> else:
> dt *= 0.8
> T_var[:] = T_var.old
> print "warning " + str(res)
>
> print datetime.now() - startTime
>
> if __name__ == '__main__':
> raw_input("Press to proceed...")
>
>
> ___
> 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
Re: Question regarding Boundary Condition Implementation
Daniel - Two things going on here: - Your constraint is being applied at the tip of the wedge of the cylindrical coordinates. Since this "face" has no area, that Dirichlet condition doesn't generate any flux in or out of the domain. Offsetting the cylinder slightly from the origin resolves this: mesh = CylindricalGrid1D(Lr = L,nr=nx) + [[L/nx]] I recommend giving some thought as to whether it even makes sense to apply a Dirichlet condition along the centerline of a cylinder. - Because you've introduced coefficients which change with temperature, your equation is now nonlinear and so sweeping is required (2 is apparently enough). At each time step, execute: for sweep in range(2): res = eq.sweep(dt=dt) rho.setValue(rho_l,where = PCM & (T>350) ) k.setValue(k_l,where = PCM & (T>350)) Cp.setValue(Cp_l,where = PCM & (T>350)) With these changes, the evolution seems generally similar to your Grid1D case. - Jon > On Apr 15, 2019, at 1:28 PM, Daniel DeSantis wrote: > > Jonathan, > > I have just a few follow up questions. Namely, regarding to working within > cylindrical coordinates. If I switch the mesh to cylindrical, it seems the > temperature profile never changes, unless I go to extremely small time steps. > That strikes me as odd. What do you think? > > Also, the shifting profile seems to reverse, with the left boundary condition > dropping to ~300K over time while in Cartesian coordinates, the temperature > rises to 600K over time (a result I would expect). I've included a copy of > the work with the changes you've suggested previously. If you look at lines > 23/24, you can toggle the cylindrical grid on and off, and in lines 30-34, > I've included options for the time steps that work for each respective > coordinate system. Do you have any suggestions? > > Finally, in reference to the velocity I would like to incorporate, I > understand the difference between Cell Variables and Face Variables, but why > is the rank of the velocity -1? What does the Face Variable rank indicate? > > Thank you, > Dan > > On Thu, Apr 11, 2019 at 8:14 AM Guyer, Jonathan E. Dr. (Fed) via fipy > wrote: > > > > On Apr 10, 2019, at 4:42 PM, Daniel DeSantis wrote: > > > > 1) Set a different Cp, rho, and k in just the PCM domain when the > > temperature in that domain crosses over a certain melt temperature (350K). > > I tried an if statement inside the solving loop and I think that has > > worked, based on some print testing I did in the loop. I just wanted to get > > your opinion and see if that was the right way to go. > > Seems reasonable, as long as you're calling `.setValue()` and not redefining > Cp, rho, and k. > > In fact, I wouldn't use an `if`: > > >>> for step in range(steps): > ... T.updateOld() > ... eq.solve(dt=dt) > ... rho.setValue(rho_melt, where=PCM & (T > 350)) > > > > 2) I'd like to move this to cylindrical coordinates with the lengths being > > radii. I'm not sure if the cylindrical grid is working though. And if it > > isn't, is there a way to convert the diffusion equation into something that > > would solve the problem in Cartesian? More specifically, how would you > > write an appropriate equation for the heat transfer of a radial system in > > Cartesian coordinates? In short, how would you write this: > > > > in FiPy? > > The CylindricalGrids work fine, except for taking the gradient of the field, > which is an unusual need. > > > > > 3) Ideally, I would move this to a 2D system which is the most confusing to > > me. In a 2D system, heat would be transmitted radially while air would be > > flowing axially. The air would cool as it passes through the tube. The > > diffusion term in the axial direction would be significantly lower than the > > convection term. Can I use the same heat transfer equation you've suggested > > slightly modified? I would guess the equation to be something like this: > > > > eq = TransientTerm(coeff = rho*Cp,var = T) + > > PowerLawConvectionTerm(coeff=rho*Cp*v,var = T) == DiffusionTerm(coeff=k,var > > = T) > > > > s.t. v = velocity of the air. > > Looks OK > > > I would also think that I would have to set the value of v to 0 at the wall > > and beyond, which means v would be a CellVariable like rho and k. > > Actually, v should be a rank-1 FaceVariable. > > Anisotropic diffusion is supported in FiPy (see > https://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.anisotropy.html), > but are you really sure you have lower axial diffusivity? It seems to me > more likely that transport is simply dominated by axial flow. > > > I also guess this equation would be subject to any changes that would be > > made in question 2, regarding writing a cylindrical equation in a cartesian > > system. What do you think? > > Use a CylindricalGrid and don't change the equation at all. > > > Again, thank you so much for your help and I hope I
