Re: Some questions on the viewer
On Thu, Jan 23, 2020 at 11:37 AM Guyer, Jonathan E. Dr. (Fed) via fipy wrote: > > I'm not sure why you are skeptical. FiPy does not use VTK internally, so the > fact that the FiPy mesh is defined appropriately for its own use has nothing > to do with whether we export it correctly as VTK. It would appear that we do > not output VTK_CONVEX_POINT_SETs correctly. FiPy doesn't order the cell to face indices, face to vertex indices or cell to vertex indices in any particular way. It doesn't matter for the FV method. It only matters for the direction of the normals and that's accounted for. I remember that those orderings are wrong for VTK as that does require a specific ordering (right (or left) hand rule or something like that). -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: CylindricalGrid1D mesh volumes
On Fri, Jan 17, 2020 at 5:31 AM Pavel Aleynikov wrote: > > sort of ratio so the absolute value of the volume won't matter. > Unless you have a source with a defined integral absolute value. But I agree > that in general "theta = 1" is not a problem, as long as it's clear from the > documentation. The documentation is probably not clear on this, but thanks for raising the issue. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: CylindricalGrid1D mesh volumes
On Wed, Jan 15, 2020 at 4:21 AM Pavel Aleynikov wrote: > > Hi, > > How are "mesh.cellVolumes" defined in fp.CylindricalGrid1D case? > The surface of a circle grid is not equal pi*r^2. > > import fipy as fp > L= 1.; nx = 1000 > mesh = fp.CylindricalGrid1D(dx=L/nx, nx=nx) > print(mesh.cellVolumes.sum()) > >> 0.5 > > Why not pi? I'm not sure. It's an arbitrary choice though. The angle was chosen as "1" rather than "2 * pi". The volume of an element is "theta * r * dr" where "theta" is the angle, "r" are the cell centers and "dr" is the cell spacing. It's possible that by choosing "theta=1", then the "theta" can be omitted saving an extra operation. > How should I integrate Variables on such a grid? > (2*pi*Var*mesh.cellVolumes).sum()? Makes sense. Argulably, these quantities should be calculated in some sort of ratio so the absolute value of the volume won't matter. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Query on AdsorbingSurfactantEquation.py consumptionCoeff
Feel free to add any relevant comments to the issue, https://github.com/usnistgov/fipy/issues/690. On Mon, Jan 13, 2020 at 1:36 PM Daniel Wheeler wrote: > > Hi Chaitanya, > > I think you're correct. It seems that term should be multiplied by the > time step. The tests all pass with your change. Would you like to > submit a pull-request for that change? > > The way time stepping works is a bit weird in that equation so it > wasn't clear to me what my intention was when I wrote that code. > Anyway, I'll add it as an issue in FiPy. > > Thanks! > > On Thu, Jan 9, 2020 at 9:13 PM Chaitanya Joshi > wrote: > > > > Hi, > > > > I am working with the leveler.py script in the level set examples. The > > script chooses dt value for solving the equations based on the maximum > > velocity of the interface at any given instant and the CFL number. > > > > Script: leveler.py > > dt = cflNumber * cellSize / extOnInt.max() > > > > I changed that to a user specified fixed dt. To test the changes I made, I > > ran the script for a flat interface and plotted acceleratorVar.interfaceVar > > as a function of time for different dt values. Unfortunately, I was unable > > to get dt convergence. By trial and error, I found out that if I call > > AdsorbingSurfactantEquation with consumptionCoeff=None, > > I get dt convergence. I then went through the > > AdsorbingSurfactantEquation.py to see if I can find something. I think > > there is a dt term missing in the AdsorbingSurfactantEquation.py. > > > > Script: AdsorbingSurfactantEquation.py > > self.eq = TransientTerm(coeff = 1) - > > ExplicitUpwindConvectionTerm(SurfactantConvectionVariable(distanceVar)) > > self.dt = Variable(0.) > > mesh = distanceVar.mesh > > adsorptionCoeff = self.dt * bulkVar * rateConstant #dt > > multiplication > >spCoeff = adsorptionCoeff * distanceVar._cellInterfaceFlag > > scCoeff = adsorptionCoeff * distanceVar.cellInterfaceAreas / > > mesh.cellVolumes > > self.eq += ImplicitSourceTerm(spCoeff) - scCoeff > > if otherVar is not None: > > otherSpCoeff = self.dt * otherBulkVar * otherRateConstant * > > distanceVar._cellInterfaceFlag #dt multiplication > > otherScCoeff = -otherVar.interfaceVar * scCoeff > > self.eq += ImplicitSourceTerm(otherSpCoeff) - otherScCoeff > > > > if consumptionCoeff is not None: > > self.eq += ImplicitSourceTerm(consumptionCoeff) #This is > > where I think a dt term is missing > > > > I believe the line should be > > if consumptionCoeff is not None: > > self.eq += ImplicitSourceTerm(consumptionCoeff*self.dt) > > > > Please let me know if I am thinking on the right lines. > > > > Regards, > > Chaitanya > > > > > > > > ___ > > fipy mailing list > > fipy@nist.gov > > http://www.ctcms.nist.gov/fipy > > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > -- > Daniel Wheeler -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Query on AdsorbingSurfactantEquation.py consumptionCoeff
Hi Chaitanya, I think you're correct. It seems that term should be multiplied by the time step. The tests all pass with your change. Would you like to submit a pull-request for that change? The way time stepping works is a bit weird in that equation so it wasn't clear to me what my intention was when I wrote that code. Anyway, I'll add it as an issue in FiPy. Thanks! On Thu, Jan 9, 2020 at 9:13 PM Chaitanya Joshi wrote: > > Hi, > > I am working with the leveler.py script in the level set examples. The script > chooses dt value for solving the equations based on the maximum velocity of > the interface at any given instant and the CFL number. > > Script: leveler.py > dt = cflNumber * cellSize / extOnInt.max() > > I changed that to a user specified fixed dt. To test the changes I made, I > ran the script for a flat interface and plotted acceleratorVar.interfaceVar > as a function of time for different dt values. Unfortunately, I was unable to > get dt convergence. By trial and error, I found out that if I call > AdsorbingSurfactantEquation with consumptionCoeff=None, > I get dt convergence. I then went through the AdsorbingSurfactantEquation.py > to see if I can find something. I think there is a dt term missing in the > AdsorbingSurfactantEquation.py. > > Script: AdsorbingSurfactantEquation.py > self.eq = TransientTerm(coeff = 1) - > ExplicitUpwindConvectionTerm(SurfactantConvectionVariable(distanceVar)) > self.dt = Variable(0.) > mesh = distanceVar.mesh > adsorptionCoeff = self.dt * bulkVar * rateConstant #dt multiplication >spCoeff = adsorptionCoeff * distanceVar._cellInterfaceFlag > scCoeff = adsorptionCoeff * distanceVar.cellInterfaceAreas / > mesh.cellVolumes > self.eq += ImplicitSourceTerm(spCoeff) - scCoeff > if otherVar is not None: > otherSpCoeff = self.dt * otherBulkVar * otherRateConstant * > distanceVar._cellInterfaceFlag #dt multiplication > otherScCoeff = -otherVar.interfaceVar * scCoeff > self.eq += ImplicitSourceTerm(otherSpCoeff) - otherScCoeff > > if consumptionCoeff is not None: > self.eq += ImplicitSourceTerm(consumptionCoeff) #This is where > I think a dt term is missing > > I believe the line should be > if consumptionCoeff is not None: > self.eq += ImplicitSourceTerm(consumptionCoeff*self.dt) > > Please let me know if I am thinking on the right lines. > > Regards, > Chaitanya > > > > ___ > 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 ]
Re: Re[2]: Need help to set up the model
Hi Mohammad, I can't actually run the example because I don't have netCDF4 installed. Could you possible formulate a simpler example without that requirement? There are a number of things that I would change in the example, but I need to be able to run it to see if those are relevant to the issues that you're seeing. Cheers, Daniel On Sun, Nov 24, 2019 at 4:02 PM Mohammad Madani wrote: > > Thank you for helping me with the mesh. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Fluid Dynamic Basics
Hi Daniel,
Solving the flow equations isn't easy and we certainly don't provide
much guidance on how to solve them with FiPy. We do have the Stokes
cavity example [1] and reactive wetting example [2], which might be a
good place for you to start. The Stokes cavity example demonstrates
the simple algorithm [3] which basically turns the continuity equation
into an equation that solves for pressure. However, this becomes more
complicated if the convection terms are included in the momentum
equation.
I have solved compressible flow problems with FiPy for my research.
FWIW, I have code for this, see [4, 5]. However, I don't think that
the code is much use to you other than as a reference for what is
possible with FiPy given enough persistence.
You might also want to try Dolfyn [6]. It might have an example that
you can use right out of the box.
Cheers,
Daniel
[1]:
https://www.ctcms.nist.gov/fipy/examples/flow/generated/examples.flow.stokesCavity.html
[2]:
https://www.ctcms.nist.gov/fipy/examples/reactiveWetting/generated/examples.reactiveWetting.liquidVapor1D.html#module-examples.reactiveWetting.liquidVapor1D
[3]: https://en.wikipedia.org/wiki/SIMPLE_algorithm
[4]: https://gist.github.com/wd15/c28ab796cb3d9781482b01fb67a7ec2d
[5]: https://journals.aps.org/pre/abstract/10.1103/PhysRevE.82.051601
[6]: https://www.dolfyn.net/index_en.html
On Thu, Aug 29, 2019 at 10:33 AM Daniel DeSantis wrote:
>
> Hello to the FiPy Team:
>
> I was recently reviewing some old work and thinking of converting it over to
> FiPy and realized that, while I had done some successful work in FiPy (with
> your assistance several times), there are a few gaps in my knowledge of how
> to apply FiPy to some CFD type problems. I apologize for asking questions you
> may feel are clear in your guide, but I am trying to learn how to use your
> awesome program a bit better so I don't have to keep asking questions.
>
> I'm essentially looking at coupling the continuity equation for a fluid with
> the equations for motion or the equations of change. So, we start with this:
>
> \nabla(\rho*v) = -\frac{\partial \rho}{\partial t}
>
> But frequently I have a steady state situation in which:
> \nabla(\rho*v) = 0
>
> In either case, I'm not sure how the continuity equation up in FiPy. Do you
> have any suggestions? Should I set up multiple convection terms like this:
>
> PowerLawConvectionTerm(coeff=rho, var = v_x) +
> PowerLawConvectionTerm(coeff=rho, var = v_y) +
> PowerLawConvectionTerm(coeff=rho, var = v_z) = 0
>
> Similarly, we have momentum equations as:
>
> \nabla(\rho v_x) = -\frac{\partial{P}}{\partial x} +\rho g_x + \mu \nabla^2v_x
>
> \nabla(\rho v_y) = -\frac{\partial{P}}{\partial y} +\rho g_y + \mu \nabla^2v_y
>
> \nabla(\rho v_z) = -\frac{\partial{P}}{\partial z} +\rho g_z + \mu \nabla^2v_z
>
> I frequently am able to ignore the \rho g_i term but if I didn't, I believe I
> could treat this as a source term in Fipy, and simply add it to my equation.
> I understand how to do the diffusion and convection terms, from your guides.
>
> However, I don't know what to do with the pressure term. I don't usually have
> a consistent pressure in all dimensions so I would think that the pressure
> term should be a cell variable, but I don't know how to add it. I tried
> treating it as P.faceGrad, but I get errors. I've tried it as a
> PowerLawConvectionTerm, and I don't seem to get a result. Perhaps the problem
> is with the continuity equation?
>
> As a sort of classic example, I've uploaded a code that describes two plates
> with a fluid in between. In it, fluid is flowing between two plates. The
> bottom plate is stationary, the top one is moving. There's a pressure
> difference over the length of the plates. It's at steady state, but I'll
> probably end up playing with it as a non-steady state system eventually. I've
> tried it before, and perhaps I'm graphing it incorrectly, but it seems like
> it just presents a steady state solution.
>
> If I've got this completely wrong, could you suggest how to set up a
> continuity equation coupled with a momentum equation and, perhaps how to
> handle a first order differential term like the pressure gradient?
>
> Thank you,
>
> --
> Daniel DeSantis
>
>
> _______
> 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 ]
Re: Solution to the sine-gordon equation
On Thu, Mar 21, 2019 at 1:14 PM Meier Quintin wrote: > > Dear Daniel, > Thanks a lot for the help. I got it to work by adding a transientterm(). > > If I try to use your second trick though: > > eq = (DiffusionTerm(coeff=(1.0), var=phi) - numerix.sin(phi) + phi * > numerix.cos(phi)-ImplicitSourceTerm(numerix.cos(phi)) == > TransientTerm(var=phi)) I think this has the parentheses in the wrong place. Does the code that I included in my last email work for you as is? This (without the line break) eq = TransientTerm() == DiffusionTerm() - numerix.sin(phi) + phi * numerix.cos(phi) - ImplicitSourceTerm(numerix.cos(phi)) > I get the error: > > ExplicitVariableError: Terms with explicit Variables cannot mix with Terms > with implicit Variables. That's because of the parentheses I think. > Also, I dont fully understand the expression: I assume the implicit term is > to neutralize the phi * numerix.cos(phi), so should it not be: > phi*ImplicitSourceTerm(numerix.cos(phi))? ImplicitSourceTerm includes the extra variable multiplier automatically. It put the source into the matrix diagonal to stabilize the solution. > Also: if I would want a second derivative in time in there to solve the > time-dependent problem, would the best way to do that to define a second > variable: > > dphidt=TransientTerm(var=phi) > > and solve the coupled equations or is there a different way? Yes, right, it can be split into two hyperbolic equations, but FiPy doesn't do well with those types of equations. It doesn't have the correct numerical schemes. I think CLAWPACK might be a better bet for that. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Solution to the sine-gordon equation
On Thu, Mar 21, 2019 at 11:48 AM Daniel Wheeler wrote: > > with a transient term in the equation. For small phi, this should > never be negative, since phi should only grow, but for delta_t > 1 That's wrong, phi should shrink, but not go negative. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Solution to the sine-gordon equation
Hi Qunitin,
I think you need to relax the updates on this. Introducing a transient
term does that. I think that the b vector causes and instability. The
b vector is
phi_old - delta_t * sin(phi_old)
with a transient term in the equation. For small phi, this should
never be negative, since phi should only grow, but for delta_t > 1
this can be negative. It's the same reason that explicit diffusion has
a time step restriction.
You can get around this restriction by linearizing the source term so
that the source looks like
- numerix.sin(phi) + phi * numerix.cos(phi) -
ImplicitSourceTerm(numerix.cos(phi))
There is now no time step restriction since the b vector is just
"phi_old" when phi is small.
After 100 sweeps, the non-linearized source is at a residual of 0.0022
(with dt=1) and the linearized source is at 0.0009 (with dt=1000),
which is somewhat better.
Possibly using Newton iterations could improve the convergence rate further.
Hope that helps.
Cheers,
Daniel
import numpy as np
from fipy import *
nx = 1
dx = 0.01
mesh = Grid1D(nx=nx, dx=dx)
phi = CellVariable(mesh=mesh, name="phi", hasOld=True)
phi.constrain(2*np.pi, where=mesh.facesLeft)
phi.constrain(0., where=mesh.facesRight)
# requires dt < 1
#eq = TransientTerm() == DiffusionTerm() - numerix.sin(phi)
# no dt requirment
eq = TransientTerm() == DiffusionTerm() - numerix.sin(phi) + phi *
numerix.cos(phi) - ImplicitSourceTerm(numerix.cos(phi))
sweeps=1000
view = Viewer(phi)
for i in range(sweeps):
res = eq.sweep(var=phi, dt=1000)
print(res)
view.plot()
phi.updateOld()
raw_input('stopped')
On Thu, Mar 21, 2019 at 7:26 AM Meier Quintin wrote:
>
> Dear fipy community,
> I’m new to FiPy, so forgive me if this is a trivial question.
> I’m trying to find the kink solution to the (static) sine-gordon equation
> using this simply piece of code.
>
> import numpy as np
> from fipy import *
> nx = 1
> dx = 0.01
> mesh = Grid1D(nx=nx,dx=dx)
> phi = CellVariable(mesh=mesh, name=“Phi")
> phi.constrain(2*np.pi, where=mesh.facesLeft)
> phi.constrain(0., where=mesh.facesRight)
> eq = (DiffusionTerm(coeff=(1.0), var=phi) - numerix.sin(phi) == 0.0)
> sweeps=10
> for i in range(sweeps):
> eq.sweep(var=phi)
> MatplotlibViewer(phi)
>
> Unfortunately, I do not manage to get a numerically stable solution, the
> first sweep gives me a straight line and additional sweeps lead to numerical
> instabilities without convergence.
> I would be very thankful if someone could hint me to what I’m doing wrong.
>
> Best wishes,
> Quintin Meier
>
> ___
> 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 ]
Re: is it easy to get fipy to make full use of dual CPUs in a single workstation?
Hi Drew, My advice would be to try FiPy in parallel on an existing device and do some speed tests to ensure you're happy. However, there's nothing about a dual CPU workstation that should be an issue. To use FiPy in parallel, you will need to have Trilinos installed. On Sun, Feb 10, 2019 at 4:24 PM Drew Davidson wrote: > > 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. See below. > 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. That doesn't mean that FiPy is being efficient. Most likely, the solvers are running in threaded mode and FiPy is not really running in parallel. To ensure that FiPy is actually running in parallel, try running the test in parallel, https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#parallel-tests. Also try running this script, https://github.com/usnistgov/fipy/blob/develop/examples/parallel.py using mpirun -np 4 python parallel.p to ensure that you have Trilinos installed correctly. Finally, once you have those working, then try testing your specific example using PySparse on a single node and then using Trilinos on multiple cores to ensure that you're getting a speed up. When running with a single core, you should only see one core being used. If all are being used then you may have threading turned on. This needs to be switched off when compiling Trilinos. Cheers, Daniel -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Boussinesq Equations
It certainly won't be that simple. You'll need to look through the
literature to figure it out. I suspect it will require re-introducing
a pressure into the equations and then transforming the continuity
equation into an equation for pressure.
On Thu, Oct 11, 2018 at 4:06 AM fgendr01 wrote:
>
> Thanks Daniel for your answer,
> I used the continuity equation (nabla vector velocity) to obtain the equation
> in temperature.
> But actually, maybe we should use it as a constraint to my problem.
> But now, how can I create this constraint ? with a new equation ?
> I tried : velocity.divergence == 0 but it doesn’t work.
>
> Thank you,
>
> Fabien
>
>
>
>
> Le 10 oct. 2018 à 17:43, Daniel Wheeler a écrit :
>
> Don't you still have a $\nabla . \vec{u} = 0$ equation though? It
> doesn't go away. That equation becomes like a constraint.
>
> https://www.comsol.com/multiphysics/boussinesq-approximation
>
> On Wed, Oct 10, 2018 at 5:58 AM fgendr01 wrote:
>
>
> Hi Daniel,
> Thank you for your answer.
> I thank you for trying to solve my problem.
> About my set of the equation here is my reasoning.
>
>
>
> --
> 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 ]
--
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: cylindrical coordinates
On Thu, Sep 20, 2018 at 4:15 AM Martinus WERTS wrote: > > Now I would like to go 3D, and the symmetry of our system would allow to > use a 2D cylindrical grid (r,z) - with zero flux at z=0 and z=L (and > r=0), and either Dirichlet/zero flux at r=R. Looking at the mailing list > archive and GitHub, it appears that cylindrical coordinates are at the > moment not working properly (missing factor) Yes, that's correct. Sorry about that. > Is this still the case? > > If so, I will start with a simple 3D Cartesian mesh, and then perhaps > move to more adapted meshes, depending on the calculation speed and the > precision required. Perhaps the cylindrical coordinates will be fixed in > a near future? Maybe, the developers of FiPy aren't doing an awful lot of work on it so there is no timeline for fixing it. > Many thanks for any advice that you may have Sorry that I can't be more helpful. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Boussinesq Equations
Don't you still have a $\nabla . \vec{u} = 0$ equation though? It
doesn't go away. That equation becomes like a constraint.
https://www.comsol.com/multiphysics/boussinesq-approximation
On Wed, Oct 10, 2018 at 5:58 AM fgendr01 wrote:
>
> Hi Daniel,
> Thank you for your answer.
> I thank you for trying to solve my problem.
> About my set of the equation here is my reasoning.
>
--
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: interpolation taking too long in an adaptive meshing effort (remeshing)
On Sun, Aug 12, 2018 at 2:52 PM Drew Davidson wrote:
>
> 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?
The call method to cell variable takes a nearestCellIDs argument which
can greatly speed up the interpolation if you're interpolating between
the same meshes multiple times. For adaptive meshing the meshes need
to retain knowledge of the adaptation to make the interpolation faster
otherwise a brute force technique is required. Scipy or Sklearn
probably has better algorithms for finding nearest neighbors so that
could be used to find the nearest neighbors instead of FiPy's method.
To use the nearestCellIDs argument, try this,
import fipy
m0 = fipy.Grid2D(nx=10, ny=10)
m1 = fipy.Grid2D(nx=20, ny=10)
v0 = fipy.CellVariable(mesh=m0, value= m0.x * m0.y)
v1 = fipy.CellVariable(mesh=m1)
# very slow and can be replaced
ids = m0._getNearestCellID(m1.cellCenters)
# should be fast at nearest neighbors are already calculated
v1[:] = v0(points=m1.cellCenters, order=1, nearestCellIDs=ids)
print v1
You could replace "_getNearestCellID" sklearn's,
http://scikit-learn.org/stable/modules/neighbors.html, without doing
anything to fipy. "ids" is just a numpy array.
> 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?
To get data out of fipy, you don't need to do much more than,
import pandas
pandas.DataFrame(dict(x=m1.x, y=m1.y, v=v1)).to_csv('test.csv')
for example. This will write the position and values of a cell
variable to a CSV file.
--
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy Heat Transfer Solution
Apologies for the slow reply. I hadn't noticed your response until now. See below. On Wed, Aug 8, 2018 at 1:09 PM Daniel DeSantis wrote: > > Daniel, > > Thank you very much for your help.The code does run much better. > > I was curious about these lines of code. I'm hoping to understand what this > means so that I can use it if needed, next time. > > Thank you again, > Dan DeSantis > > constraint_value = FaceVariable(mesh=mesh) > T.faceGrad.constrain(constraint_value,mesh.facesLeft) > > for sweep in range(sweeps): > constraint_value[:] = h * ((60. + 273.15) - T.faceValue) I only used the above code as it happened to work. It should work without having to do the update in the loop, but something isn't right with the way FiPy is updating constraints in the case when the constraint is not a bare FaceVariable. In other words, we make a face variable and then say that the constraint depends on that face variable. We have to explicitly update the face variable in the loop. The expression "h * ((60. + 273.15) - T.faceValue" should also be a face variable and so we should be able to use "T.faceGrad.constrain(h * ((60. + 273.15) - T.faceValue, mesh.facesLeft)" and not update in the loop. That doesn't work, however. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy Heat Transfer Solution
Hi Daniel, I updated your code a bit and it seems to do better. See https://pastebin.com/51xfqWH3 for the full version. I changed the "eqX" equation to param = k_c * (P - Peq) maxX = 0.9 largeValue = 1e9 constraint = largeValue * (X > maxX) big_value = 1e9 eqX = TransientTerm(var = X) == param - ImplicitSourceTerm(param, var=X) + maxX * constraint - ImplicitSourceTerm(constraint, var=X) to have a constraint on the maximum value of X. I'm not sure whether this is physical, but it does constrain the max value of X. I changed the "eqY" equation to be, eqT = TransientTerm(rho_b_MH * Cp_b, var=T) == DiffusionTerm(coeff = L_eff, var=T) + rho_b_MH*(dH/M)*rc to have the transient coefficient inside the term. Again, not sure if that matters, but did it anyway. I changed the constraint to be constraint_value = FaceVariable(mesh=mesh) T.faceGrad.constrain(constraint_value,mesh.facesLeft) and then updated the constraint in the loop with for sweep in range(sweeps): constraint_value[:] = h * ((60. + 273.15) - T.faceValue) If FiPy worked correctly that wouldn't be necessary, but evidently it doesn't. Anyway that all seems to work at least as far as I can tell. Cheers, Daniel -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy Heat Transfer Solution
On Wed, Jul 18, 2018 at 4:36 PM, Daniel DeSantis wrote: > > > If I wanted to run this in a cylindrical coordinate grid, the equation would > become: > > > Can I just call the following in FiPy to account for this coordinate system > change?: > > mesh = CylindricalGrid1D(nx = 100.,Lx = 1., dx = 1./100.) Yes, that should work. There is a bug (see https://github.com/usnistgov/fipy/issues/547) when calculating gradients, but I don't think that applies to diffusion terms. Jon, is that correct? > If not, how do I account for the dependent variable when it is not in the > differential? Essentially, how do can I put lambda/r in the coefficient etc.? You can multiply through by r and then split the transient term into an implicit TransientTerm and a source term. It gives an extra source term that is dependent on the time step. > 2) Is the correct way to add a source term simply to write something like: > > eq = TransientTerm() == DiffusionTerm() + S > > or > > eq = TransientTerm() == DiffusionTerm() + ImplicitSourceTerm(S)? Either one. ImplicitSourceTerm implies that the source term is "var * S", not "S" only, but the "var" is the variable that you're solving for. To use ImplicitSourceTerm the S should be negative. In your case, the source does not depend on temperature so you can't use an ImplicitSourceTerm. > 3) Regarding boundary conditions, I am trying to use a boundary condition > that has my dependent variable (T) in the boundary condition. The condition > is dT/dr = h(T - Tc), where Tc is a constant. Would this be the correct way > to do that? I'm mostly inquiring about writing the dependent variable as > T.faceValue... > > T.faceGrad.constrain(h*(T.faceValue - Tc), mesh.facesLeft) I think that should work. Test it carefully. It's never obvious to me what the sign should be. Don't trust it though. -- Daniel Wheeler ___ 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
On Mon, Jul 23, 2018 at 1:06 PM, Drew Davidson wrote: > > 'dx' is actually the thickness of the thermal contact region. It is not the > volume of the cell. If it were a volume, the ratio dx/mesh._cellDistances > would not be dimensionless. I was seeking a thermal contact region of zero > thickness (dx=0), giving a true thermal contact resistance. Sorry, my mistake on that one. > I let this thread drop because the Salazar paper I referenced gave me > expressions for effective thermal conductivity and effective thermal > diffusivity, and this allowed me to arrive at a reasonable-looking FiPy > result for transient thermal behavior. My result still deviates a bit from > the exact solution at steady state, and this deviation persists despite mesh > refinement, but I am busy with other things and have called it good enough > for now. Thanks for following up and sharing your code. -- Daniel Wheeler ___ 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
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 ]
Re: Integral term
On Thu, Jul 12, 2018 at 9:59 AM, Pavel Aleynikov wrote: > How should such equation be arranged? fp.numerix.cumsum seems to return a numpy array so won't automatically update when f1 changes. Probably best to update in the loop like this, import fipy as fp mesh = fp.Grid1D(dx=0.002, nx=100) mesh = mesh + [0.002] f1 = fp.CellVariable(name = "solution",mesh = mesh,value = 1.,hasOld=True) x = mesh.x M1 = lambda f1: fp.numerix.cumsum(f1*mesh.cellVolumes) M2 = lambda f1: fp.numerix.cumsum(x**2*f1*mesh.cellVolumes) Conv_pf = fp.CellVariable(mesh=mesh, rank=1) Diff_pf = fp.CellVariable(mesh=mesh) eq_steady1 = fp.DiffusionTerm(coeff=Diff_pf) + fp.PowerLawConvectionTerm(coeff=Conv_pf) f1.constrain(1,mesh.facesLeft) for sweep in range(10): Conv_pf[:] = M1(f1) / x**2 Diff_pf[:] = M2(f1) / x**3 eq_steady1.sweep(var=f1) print f1 f1.updateOld() -- Daniel Wheeler ___ 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
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 ]
Re: thermal contact between two regions
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 ]
Re: Moving boundary in FiPy
On Mon, Jul 2, 2018 at 3:16 PM, Zhekai Deng wrote: > Hi Daniel, > > Thanks for the pointed sources. I think I still need some help on how to > setup my boundary at irregular triangular gmsh. Even though the locations of > my boundary conditions are simple (two straight lines inclined with fixed > degree), my mesh generated from gmsh is triangular and irregular. Thus, my > boundary location lines are not fully aligned with the face edges in the > mesh. I wonder how to setup the "where" option in case like this? It depends if the boundary condition is defined with a sharp or diffuse formulation. If you're using a sharp boundary formulation then the mesh faces need to be aligned along the boundary. In that case I don't know any other way. Gmsh will certainly let you do that. > Right now, I am thinking that I need to mark my boundary by myself. I have > been trying a lot of statements like following but it seems doesn't work > out: > phi.faceGrad.constrain([a*local_shear_rate*(1-phi.faceValue)*phi.faceValue], > where = np.abs(yFace -(leftheight - (L + > xFace)*np.tan(25.0/180.0*np.pi))) I am not sure whether It is better to formulate my boundary condition as a > source term instead of using faceGrad.constrain... Exactly, formulate the boundary condition as a source and variable diffusion coefficient. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Moving boundary in FiPy
Hi Zhekai, I don't think there are any examples of using FiPy with a moving mesh. It would definitely not be easy to set up and implement with FiPy. It would also be extremely slow in FiPy due to the upfront cost of creating meshes. FiPy has been used with both phase field and level set examples (fixed mesh, Eulerian approach). FiPy has also been used for examples with the equations recast into a moving frame of reference, but this only works for simple geometry changes. Cheers, Daniel On Wed, May 9, 2018 at 5:30 PM, Zhekai Deng <zhekaideng2...@u.northwestern.edu> wrote: > Hi All, > > I wonder is there any way to implement moving boundary (with given known > mesh moving velocity) in FiPy? I couldn't find any example related to this > online and wonder if anyone might have tried it before. I basically would > like to solve an advection-diffusion equation in a moving mesh. > > Thanks! > > Zhekai > > ___ > 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 ]
Re: how to compute interface velocity in examples.phase.anisotropy?
Hi Drew,
The velocity vector is the velocity magnitude multiplied by the
normal. The normal is the gradient over the gradient magnitude so,
\vec{v} = \frac{\nabla \phi}{\nabla \phi \cdot \nabla \phi}
\frac{\partial \phi}{\partial t}
In Fipy you can get the gradient and gradient magnitude with the
operations ".grad" and ".grad.mag" on cell variables. See, for
example,
https://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.impingement.mesh40x1.html,
where a .grad.mag is used.
Cheers,
Daniel
On Tue, Apr 3, 2018 at 4:55 PM, Drew Davidson <davidson...@gmail.com> wrote:
> 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 <davidson...@gmail.com>
> 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 ]
>
--
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Spline interpolation and fipy variable
On Fri, Jan 12, 2018 at 12:42 PM, Clara Maurel <cmau...@mit.edu> wrote: > Hi Daniel, > > Thank you very much! The code runs for me as well: fantastic! > > I have to admit that the distinction between cell and face variables remains > a bit obscure to me… The domain is divided into cells or boxes in 2D. The faces are where the boxes meet and the cell centers are located at the center of the boxes. FiPy calculates/resolves fluxes between the boxes. > The cell variable X_var is the composition of the > system I am modelling. Then "G" in d2G is the Gibbs free energy (called f > here > https://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html#module-examples.cahnHilliard.mesh2DCoupled), > which is itself a function of the composition X_var (the spline function). > I don’t know if these info are more useful to figure out if d2G should be > cell or face variable? As I understood it, because I want to solve for > X_var, this should be a cell variable and nothing else. X_var is a cell variable but the diffusion coefficient can either be a cell variable or a face variable in FiPy. FiPy converts the diffusion coefficient to a face variable if it's defined as a cell variable in the input file. FiPy calculates the flux between cells so the diffusion coefficient is only used at the faces where the flux is calculated. d2G has values which are defined over the domain. How do those values correspond to a location in that domain? Why does d2G have a shape of (1001,) when there are only 1000 cells? > By higher dimension do you also mean 2D? I know 2D requires a lot of > changes, but I thought 2D required only a change of grid. In any case, 2D is > not my top priority at the moment. Ok, this is only being used for 1D right now so you're fine. > Thank you again very much for your help! Good luck with 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 ]
Re: Spline interpolation and fipy variable
Hi Clara, The following changes made it run for me. I turned d2G into a face variable and K into a variable that is updated in the sweep loop. I also removed redefining the equation. FiPy's designed so you only need to to create the equation once. I'm not sure what the values of d2G correspond to. Which spatial location do they correspond to? Faces or Cells. I'm assuming faces since there are 1001 of them. Also, this might need to be adjusted for higher dimensions. Cheers, Daniel $ git diff diff --git a/Spinodal_fipy_FeNi_varT_runs_DEBUG.py b/Spinodal_fipy_FeNi_varT_runs_DEBUG.py index 715e79b..a82387e 100644 --- a/Spinodal_fipy_FeNi_varT_runs_DEBUG.py +++ b/Spinodal_fipy_FeNi_varT_runs_DEBUG.py @@ -289,20 +289,21 @@ def Get_spline_Chebnodes(X,N,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_f return d2G + ## THIS IS WHERE I WOULD LIKE TO HAVE A FIPY VARIABLE TYPE FOR d2G. d2G = Get_spline_Chebnodes(Xf,100,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)*1.0e-9 - +d2G_var = FaceVariable(mesh=mesh, value=d2G) ## Initial diffusion coefficient -Diff = Mob*d2G +Diff = Mob*d2G_var ## Initial gradient energy (depending on the interfacial energy we want) K = 2.0*kappa[id_mean_spino] - +K_var = Variable(K) ## Initial equation -eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) - DiffusionTerm(coeff=(Mob, K)) +eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) - DiffusionTerm(coeff=(Mob, K_var)) plt.figure(figsize=(8,6)) @@ -390,16 +391,16 @@ while time < duration: Mob = M[step_T][id_mean_T] d2G = Get_spline_Chebnodes(Xf,100,T_spino[step_T],R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)*1.0e-9 - -Diff = Mob*d2G +d2G_var[:] = d2G K = 2.0*kappa[step_T] +K_var.setValue(K) print 'Mob, Diff, K' print Mob, Diff, K print 'MAX CONCENTRATION' print max(X_var) -eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) - DiffusionTerm(coeff=(Mob, K)) + print datetime.now() - startTime (END) On Thu, Jan 11, 2018 at 2:35 PM, Clara Maurel <cmau...@mit.edu> wrote: > Hi Daniel, > > Thank you for taking the time to look at this! My main code calls several > text files and subroutine. I attach everything that is needed to run the > code: -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Spline interpolation and fipy variable
Hi Clara,
Can you post the entire working code so I can try and debug?
Cheers,
Daniel
On Wed, Jan 10, 2018 at 4:43 PM, Clara Maurel <cmau...@mit.edu> wrote:
> Hello,
>
> I asked a similar questions several months ago and someone nicely answered
> with some ideas, which however did not solve my problem. So I am just trying
> another time!
>
> I want to model the dynamics of a system using the Cahn Hilliard equation,
> with a diffusion coefficient that depends on the cell variable
> (concentration X) and time.
>
> The diffusion coefficient (Diff) is proportional to the second derivative of
> the free energy (d2G). The latter function is not continuous and I would
> like to interpolate it with a B-spline (scipy.interpolate.UnivariateSpline).
> FYI, I first started with a polynomial interpolation: the fit was not
> satisfactory in terms of the physics I want to capture, but the fipy code
> worked well.
>
> Problems come when I want to use my spline function with the cell variable
> X: this always return an array where fipy would like a fipy variable, and I
> don’t know if there would be a way to do so. Would anyone have an idea on
> how I could do that?
> The suggestion I had before was to set Diff as a cell variable and update
> the value of Diff at each time step. However, doing so the behaviour of the
> system was not physical.
>
> Diff = CellVariable(mesh=mesh)
> Diff.setValue(d2G(Xf))
>
>
> Any other idea would be greatly appreciated! I would be happy to send my
> code if necessary.
>
> Thank you in advance,
> Clara
>
>
>
>
>
> Here is the part of my code in question:
>
> ## Mesh
> L = 200.
> nx = 1000
> dx = L/nx
>
> mesh = PeriodicGrid1D(nx=nx, dx=dx)
>
> ## Variable
> X_var = CellVariable(name=r"$X_{at}$", mesh=mesh, hasOld=True)
>
>
> ## Mean initial concentration
> X_mean = float(comp)/100.0
> X_mean_txt = str(comp)
> id_mean = Get_closest_id(X,X_mean)
> id_mean_spino = Get_closest_id(X_spino_low,X_mean)
>
>
> ## Initial temperature corresponding to X_mean on spinodal boundary
> T0 = T_spino[id_mean_spino]
>
> ## Initial thermal fluctuations
> noise = GaussianNoiseVariable(mesh=mesh, mean=X_mean,
> variance=0.1).value
> X_var[:] = noise
> X_var.updateOld()
> Xf = X_var.arithmeticFaceValue
>
> ## Initial mobility
> Mob = M[id_mean_spino][id_mean]
>
> ## This interpolate the second derivative of the free energy with a B-spline
> and return Spline(Xf) WHICH IS AN ARRAY...
> d2G =
> Get_spline_Chebnodes(Xf,100,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)
>
> ## Initial diffusion coefficient
> Diff = Mob*d2G
>
> ## Initial gradient energy (depending on the interfacial energy we want)
> K = 2.0*kappa[id_mean_spino]
>
> ## Initial equation
> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
> DiffusionTerm(coeff=(Mob, K))
>
> ## Time and integration parameters
> time = (T0-T_spino[id_mean_spino])/Cr
> duration = (T0-T_spino[-2])/Cr # in s
> dt = 3.154e7 # 1 year in s
> dt_max = 3.154e13 # 1,000,000 years in s
> total_sweeps = 4
> tolerance = 1.0e-3
> step_T = id_mean_spino
>
> solver = Solver()
>
> ## Start looping
> while time < duration:
>
>
>
> res0 = eq.sweep(X_var, dt=dt, solver=solver)
>
>
>
> next_time_T = (T0-T_spino[step_T+1])/Cr
> print "Next temperature step is at: " + str(next_time_T/Myr_in_s) + "
> Myr"
> print "Temperature: "+str(T_spino[step_T]-273)
>
>
>
> while time < next_time_T:
>
>
>
> for sweeps in range(total_sweeps):
> res = eq.sweep(X_var, dt=dt, solver=solver)
>
> if res < res0 * tolerance:
> time += dt
> dt *= 1.1
> dt = min(dt_max, dt)
> X_var.updateOld()
> Xf = X_var.arithmeticFaceValue
>
>
> else:
> dt *= 0.8
> X_var[:] = X_var.old
> Xf = X_var.arithmeticFaceValue
>
> print "Arrived at a temperature step!"
> step_T += 1
> id_mean_T = Get_closest_id(X,np.mean(X_var))
>
> Mob = M[step_T][id_mean_T]
>
>
>
> d2G =
> Get_spline_Chebnodes(Xf,100,T_spino[step_T],R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)
>
>
>
> Diff = Mob*d2G
>
>
>
> K = 2.0*kappa[step_T]
>
> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
> DiffusionTerm(coeff=(Mob, K))
>
>
>
>
>
>
>
>
>
>
> ___
> 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 ]
Re: Gradient in cylindrical coordinates
Also, note that the "leastSquaresGrad" doesn't have this issue: https://gist.github.com/wd15/fc34ccb2e57602fc6f9bea96d8160f4a#file-untitled-ipynb See, https://www.ctcms.nist.gov/fipy/fipy/generated/fipy.variables.html#fipy.variables.cellVariable.CellVariable.leastSquaresGrad On Mon, Jan 8, 2018 at 1:39 PM, Daniel Wheeler <daniel.wheel...@gmail.com> wrote: > Hi Leyton, > > I think there is an error in FiPy when calculating gradients with > cylindrical grids. FiPy assumes that the cylindrical faces don't count > towards the face contributions, but they do. It needs to be fixed. It > results in a systematic error of "phi / r". I'll try and fix this in > the near future. Sorry about that. > > For further details see, > > https://gist.github.com/wd15/fc34ccb2e57602fc6f9bea96d8160f4a#file-untitled-ipynb > > I've very roughly outlined the problem at the bottom of the notebook > and a better discretization. A quick fix right now would be to use a > 3D slice rather than the cylindrical grid. I'm not sure how much this > impacts other calculations with cylindrical grids in FiPy. > > I also put this in an issue. Please do add comments if you have ideas. > > https://github.com/usnistgov/fipy/issues/547 > > Thanks for reporting this. > > Cheers, > > Daniel > > > On Fri, Jan 5, 2018 at 7:23 PM, Munoz Leyton > <leyton.mu...@hirtenberger.com> wrote: >> Dear Fipy Team, >> >> >> >> First of all I would like to thank you for your amazing work! I love working >> with fipy. >> >> >> >> I have a problem when I calculate the gradient of my variable in cylindrical >> coordinates. You can see the code I used below this line. >> >> >> >> from fipy import * >> >> from fipy import CellVariable, FaceVariable, Grid2D, TransientTerm, \ >> >> UpwindConvectionTerm, DiffusionTerm, ImplicitSourceTerm, >> CylindricalGrid2D >> >> >> >> L = 89e-3 # m >> >> diameterChamber = 13.75e-3 >> >> radius = diameterChamber/2 >> >> dz = 1e-3 >> >> dr = 0.75e-3 >> >> nz = round(L/dz) >> >> nr = round(radius/dr) >> >> #mesh = Grid2D(dx=dx, nx=nx, dy=dy, ny=ny) >> >> mesh = CylindricalGrid2D(dz=dz, nz=nz, dr=dr, nr=nr) >> >> >> >> pressure_0 = CellVariable(mesh=mesh, name='Pressure Gas', hasOld=1) >> >> pressure_0.setValue(1e5) >> >> pressure_0.grad.value >> >> >> >> The pressure is constant all over the domain so I would expect the gradient >> to be zero. just as it is in Cartesian coordinates. But I get this: >> >> >> >> array([[ 2666.6667, 888.8889, 533., ..., >> >> 205.12820513, 177.7778, 156.8627451 ], >> >>[0., 0., 0., ..., >> >> 0., 0., 0.]]) >> >> >> >> Could you please help me? Am I doing something wrong? is it a bug? >> >> >> >> I would like also to contribute to the community with an example of >> compressible flow with euler equations in 2D Cartesian coordinates. I >> validated the results with SOD shock tube and I was quite happy with it. >> With whom should I talk to review it and possibly upload it? >> >> >> >> Kind regards, >> >> Leyton >> >> >> >> ACHTUNG neue e-Mail Adresse: vorname.nachn...@hirtenberger.com Nachrichten >> an die alte E-Mail Adresse werden vorübergehend auf die neue Adresse >> umgeleitet. >> >> ATTENTION new e-mail address: vorname.nachn...@hirtenberger.com Messages to >> the old address are temporarily redirected to the new address. >> >> >> >> P Bitte prüfen Sie der Umwelt zuliebe, ob der Ausdruck dieser Mail >> erforderlich ist. / Please consider your environmental responsibility before >> printing this email. >> >> >> >> Diese Nachricht und allfaellige angehaengte Dokumente sind vertraulich und >> nur für den/die Adressaten bestimmt. Sollten Sie nicht der beabsichtigte >> Adressat sein, ist jede Offenlegung, Weiterleitung oder sonstige Verwendung >> dieser Information nicht gestattet. In diesem Fall bitten wir, den Absender >> zu verstaendigen und die Information zu vernichten. Für Uebermittlungsfehler >> oder sonstige Irrtuemer bei der Uebermittlung besteht keine Haftung. >> >> This message and any attached files are confidential and intended solely for >> the ad
Re: Gradient in cylindrical coordinates
Hi Leyton, I think there is an error in FiPy when calculating gradients with cylindrical grids. FiPy assumes that the cylindrical faces don't count towards the face contributions, but they do. It needs to be fixed. It results in a systematic error of "phi / r". I'll try and fix this in the near future. Sorry about that. For further details see, https://gist.github.com/wd15/fc34ccb2e57602fc6f9bea96d8160f4a#file-untitled-ipynb I've very roughly outlined the problem at the bottom of the notebook and a better discretization. A quick fix right now would be to use a 3D slice rather than the cylindrical grid. I'm not sure how much this impacts other calculations with cylindrical grids in FiPy. I also put this in an issue. Please do add comments if you have ideas. https://github.com/usnistgov/fipy/issues/547 Thanks for reporting this. Cheers, Daniel On Fri, Jan 5, 2018 at 7:23 PM, Munoz Leyton <leyton.mu...@hirtenberger.com> wrote: > Dear Fipy Team, > > > > First of all I would like to thank you for your amazing work! I love working > with fipy. > > > > I have a problem when I calculate the gradient of my variable in cylindrical > coordinates. You can see the code I used below this line. > > > > from fipy import * > > from fipy import CellVariable, FaceVariable, Grid2D, TransientTerm, \ > > UpwindConvectionTerm, DiffusionTerm, ImplicitSourceTerm, > CylindricalGrid2D > > > > L = 89e-3 # m > > diameterChamber = 13.75e-3 > > radius = diameterChamber/2 > > dz = 1e-3 > > dr = 0.75e-3 > > nz = round(L/dz) > > nr = round(radius/dr) > > #mesh = Grid2D(dx=dx, nx=nx, dy=dy, ny=ny) > > mesh = CylindricalGrid2D(dz=dz, nz=nz, dr=dr, nr=nr) > > > > pressure_0 = CellVariable(mesh=mesh, name='Pressure Gas', hasOld=1) > > pressure_0.setValue(1e5) > > pressure_0.grad.value > > > > The pressure is constant all over the domain so I would expect the gradient > to be zero. just as it is in Cartesian coordinates. But I get this: > > > > array([[ 2666.6667, 888.8889, 533., ..., > > 205.12820513, 177.7778, 156.8627451 ], > >[0., 0., 0., ..., > > 0., 0., 0.]]) > > > > Could you please help me? Am I doing something wrong? is it a bug? > > > > I would like also to contribute to the community with an example of > compressible flow with euler equations in 2D Cartesian coordinates. I > validated the results with SOD shock tube and I was quite happy with it. > With whom should I talk to review it and possibly upload it? > > > > Kind regards, > > Leyton > > > > ACHTUNG neue e-Mail Adresse: vorname.nachn...@hirtenberger.com Nachrichten > an die alte E-Mail Adresse werden vorübergehend auf die neue Adresse > umgeleitet. > > ATTENTION new e-mail address: vorname.nachn...@hirtenberger.com Messages to > the old address are temporarily redirected to the new address. > > > > P Bitte prüfen Sie der Umwelt zuliebe, ob der Ausdruck dieser Mail > erforderlich ist. / Please consider your environmental responsibility before > printing this email. > > > > Diese Nachricht und allfaellige angehaengte Dokumente sind vertraulich und > nur für den/die Adressaten bestimmt. Sollten Sie nicht der beabsichtigte > Adressat sein, ist jede Offenlegung, Weiterleitung oder sonstige Verwendung > dieser Information nicht gestattet. In diesem Fall bitten wir, den Absender > zu verstaendigen und die Information zu vernichten. Für Uebermittlungsfehler > oder sonstige Irrtuemer bei der Uebermittlung besteht keine Haftung. > > This message and any attached files are confidential and intended solely for > the addressee(s). Any publication, transmission or other use of the > information by a person or entity other than the intended addressee is > prohibited. If you receive this in error please contact the sender and > delete the material. The sender does not accept liability for any errors or > missions as a result of the transmission. Please destroy this message and > notify the sender. > > > ___ > 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 ]
Re: Compute total free energy of a system
You might want to multiply by the cell volumes, "mesh.cellVolumes".
On Thu, Jan 4, 2018 at 10:22 AM, Anders Ericsson
<anders.erics...@solid.lth.se> wrote:
> Hi,
>
>
> I wonder if there is a simple way to compute the total free energy of a
> system (Phase-field modeling) in FiPy?
>
>
> That is e.g.:
>
>
> \begin{equation}
> F = \int_V f(\phi, c, T) + \frac{\epsilon_{\phi}^2}{2}\nabla \phi^2 dV
> \end{equation}
>
> I figured that it would be something in accordance with:
>
> def freeEnergyVolume(phi_, c_, T, epsSq):
> return (0.5 * epsSq * (phi_.grad.mag)**2 +
> f(phi_,c_,T)).cellVolumeAverage
>
> But cellVolumeAverage wouldn't give me the full energy of the domain if I'm
> not mistaken?
--
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: 2D anisotropic unstructured mesh and non-classical FV schemes
Thanks for your interest in FiPy On Mon, Oct 16, 2017 at 6:16 PM, F Hssn <fhs...@gmail.com> wrote: > > So my questions are: > > - Does fipy handle anisotropic unstructured meshes without any > problem? (I'm specifically looking to use bamg for 2D mesh adaptation > based on interpolation error reduction/equidistribution) What is an anisotropic mesh? FiPy does "handle" unstructured meshes in that the problem will solve, however, the accuracy decreases as the non-orthogonality and non-conjunctionality increases (as you pointed out). There has been an attempt to implement terms to correct for the non-orthogonality, see https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermCorrection.py and https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermNoCorrection.py You might want to try playing with the corrected diffusion term to see how well it works. I can't find any examples of its use though. See also - https://www.mail-archive.com/fipy@nist.gov/msg03758.html - https://www.mail-archive.com/fipy@nist.gov/msg03757.html - https://www.mail-archive.com/fipy@nist.gov/msg03765.html > - If not, does fipy provide (or provide a way to implement) MPFA > (multi-point flux approximation) or some other non-classical scheme > (like HMM) that allows anisotropic unstructured meshes (as discussed > in the Drioniou paper [1]) ? I don't think I can easily answer that without a great deal of work, but I did skim over the Drioniou paper, https://arxiv.org/pdf/1407.1567.pdf. I can't foresee all the issues with implementing the various schemes but one that does come to mind is having to use neighbor's neighbor values to compute fluxes, which I think these schemes require. It would require a lot of changes to FiPy to set that up efficiently. However, I will certainly keep this paper on my list of items to research further. Thanks for highlighting it. > - Or, does fipy allow any other way to make vertex-centered control > volume calculation that can take into account anisotropy and can make > sure there are no control volume overlaps (and as a result, the > coefficients stay positive, and monotonicity is maintained, and none > of the nice properties are violated) ? No, FiPy is definitely not set up for vertex centered. The above cell-centered approach is the way to go. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Error when importing PySparse
Hi Clara, I would seriously recommend using Conda and Python 2.7 for FiPy and installing PySparse from the "guyer" conda channel. I'm no longer sure which versions of PySparse work or don't work with FiPy. I have used this Dockerfile successfully to run FiPy, https://github.com/wd15/extremefill2D-dockerize/blob/master/Dockerfile. Specifically, this stuff will help you get up and running in a Conda environment: wget https://repo.continuum.io/miniconda/Miniconda2-4.3.11-Linux-x86_64.sh bash Miniconda2-4.3.11-Linux-x86_64.sh -b -p $ANACONDAPATH conda update conda conda install libgfortran=1.0 && conda clean --all conda install matplotlib && conda clean --all conda install --channel guyer scipy gmsh && conda clean --all conda install --channel guyer pysparse openmpi mpi4py && conda clean --all conda install --channel guyer trilinos && conda clean --all ~~~ There are a lot of additional packages in the Dockerfile so you might not want to use it as is, however, it might give you some clues about what's necessary to get FiPy working even if you don't use Docker. I hope this helps. Cheers, Daniel On Thu, Jul 20, 2017 at 9:38 AM, Clara Maurel <cmau...@mit.edu> wrote: > Hello, > > I have trouble trying to install PySparse on a Linux cluster. The pip > command: > > pip install --user pysparse > > did not work. > > I found on a forum this command: > > pip install --user git+http://git.code.sf.net/p/pysparse/git#egg=PySparse > > after which it said the installation was successful. However, when I try to > import Pysparse, I have an error message that says: > > Traceback (most recent call last): > File "", line 1, in > File > "/home/cmaurel/.local/lib/python2.6/site-packages/pysparse/__init__.py", > line 13, in > from sparse import spmatrix > File > "/home/cmaurel/.local/lib/python2.6/site-packages/pysparse/sparse/__init__.py", > line 6, in > from pysparseMatrix import * > File > "/home/cmaurel/.local/lib/python2.6/site-packages/pysparse/sparse/pysparseMatrix.py", > line 57, in > from pysparse.sparse import spmatrix > ImportError: > /home/cmaurel/.local/lib/python2.6/site-packages/pysparse/sparse/spmatrix.so: > undefined symbol: __intel_security_cookie > > > > I also independently tried to run: > > python setup.py build > > but this was also unsuccessful with the error: > > creating build/temp.linux-x86_64-2.6 > creating build/temp.linux-x86_64-2.6/src > compile options: '-DLENFUNC_OK=1 -DATLAS_INFO="\"3.8.4\"" -Isrc > -I/usr/include -I/usr/lib64/python2.6/site-packages/numpy/core/include > -I/usr/include/python2.6 -c' > icc: src/spmatrixmodule.c > icc: command line warning #10006: ignoring unknown option '-fwrapv' > icc: command line warning #10006: ignoring unknown option '-fwrapv' > icc: error #10236: File not found: 'src/spmatrixmodule.c' > icc: command line error: no files specified; for help type "icc -help" > icc: command line warning #10006: ignoring unknown option '-fwrapv' > icc: command line warning #10006: ignoring unknown option '-fwrapv' > icc: error #10236: File not found: 'src/spmatrixmodule.c' > icc: command line error: no files specified; for help type "icc -help" > error: Command "icc -fno-strict-aliasing -O2 -g -pipe -Wall > -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector > --param=ssp-buffer-size=4 -m64 -mtune=generic -D_GNU_SOURCE -fPIC -fwrapv > -DNDEBUG -O2 -g -pipe -Wall -Wp,-D_FORTIFY_SOURCE=2 -fexceptions > -fstack-protector --param=ssp-buffer-size=4 -m64 -mtune=generic > -D_GNU_SOURCE -fPIC -fwrapv -fPIC -DLENFUNC_OK=1 -DATLAS_INFO="\"3.8.4\"" > -Isrc -I/usr/include -I/usr/lib64/python2.6/site-packages/numpy/core/include > -I/usr/include/python2.6 -c src/spmatrixmodule.c -o > build/temp.linux-x86_64-2.6/src/spmatrixmodule.o" failed with exit status 1 > > Would anyone have an idea on how to fix this? > > Any help would be very much appreciated! Thank you in advance! > > Clara > > > ___ > 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 ]
Re: Time varying boundary conditions
Hi Adrian, It's a little bit of a hack, but I got something to work, see below. Unfortunately, I think you need to explicitly delete the fixed value constraint from the CellVariable's faceConstraints list, which requires some knowledge of FiPy's internals and remember the order in which the constraints were added. I would be nice if "var.constrain" returned some sort of constraint object which could be deleted. Cheers, Daniel import fipy as fp mesh = fp.Grid1D(nx=10) var = fp.CellVariable(mesh=mesh, value=0.) var.constrain(1., mesh.facesLeft) var.constrain(0., mesh.facesRight) eqn = fp.TransientTerm() == fp.DiffusionTerm() viewer = fp.Viewer(var) dt = 0.1 for _ in range(100): eqn.solve(var, dt=dt) viewer.plot() del var.faceConstraints[0] var.faceGrad.constrain(1., mesh.facesLeft) for _ in range(1000): eqn.solve(var, dt=dt) viewer.plot() On Wed, Jun 21, 2017 at 2:56 PM, Adrian Jacobo <ajac...@mail.rockefeller.edu> wrote: > Hi All, > > I have a problem where I would want to specify fixed value boundary > conditions for times t conditions for t>=t0. The way I've implemented this is: > > > a.constrain(a_0,(mesh.exteriorFaces)) > > while t if t > = t0: > a.faceGrad.constrain(0.,mesh.exteriorFaces) > eq.solve(dt=timeStepDuration) > t+=dt > > But what it seems to be happening is that the fixed value bc is not > being overwritten by the fixed gradient one. How can I do to remove the > fixed value bc? Is there a way to check the boundary conditions on a > variable? > > Thanks, > Adrian. > > ___ > 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 ]
Re: Conda installation not working
Hi Adrian, Sorry about the installation issues. I've started using FiPy in a Docker instance where I fix all the version numbers of the dependencies as I've had many issues like this. See, https://github.com/wd15/extremefill2D-dockerize/blob/master/Dockerfile Some of this Dockerfile has extra packages that aren't required for FiPy, but it will give you a working version of FiPy. Cheers, Daniel On Thu, Jun 8, 2017 at 10:52 AM, Adrian Jacobo <ajac...@mail.rockefeller.edu> wrote: > Hi, > > For some reason my fipy conda installation stopped running. I've > removed it and installed it again but I keep getting the same error: > > >>> import fipy > Traceback (most recent call last): >File "", line 1, in >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/fipy/__init__.py", > line 45, in > from fipy.boundaryConditions import * >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/fipy/boundaryConditions/__init__.py", > line 2, in > from fipy.boundaryConditions.fixedFlux import * >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/fipy/boundaryConditions/fixedFlux.py", > line 38, in > from fipy.tools import numerix >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/fipy/tools/__init__.py", > line 85, in > serial, parallel = serialComm, parallelComm = _getComms() >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/fipy/tools/__init__.py", > line 77, in _getComms > from fipy.tools.comms.dummyComm import DummyComm >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/fipy/tools/comms/dummyComm.py", > line 37, in > from fipy.tools import numerix >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/fipy/tools/numerix.py", > line 76, in > import numpy as NUMERIX >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/numpy/__init__.py", > line 112, in > from ._globals import ModuleDeprecationWarning, > VisibleDeprecationWarning > ImportError: No module named _globals > > I could trace the error to Numpy: > > import numpy > Traceback (most recent call last): >File "", line 1, in >File > "/Users/ajacobo/anaconda/envs/fipyenv/lib/python2.7/site-packages/numpy/__init__.py", > line 112, in > from ._globals import ModuleDeprecationWarning, > VisibleDeprecationWarning > ImportError: No module named _globals > > As far as I could figure out anaconda updated python to Python 2.7.13 > and the numpy version that ships with fipy is incompatible with this > version. Could you update numpy in the anaconda installation? I tried to > force numpy to update but couldn't. > > Thanks! > Adrian. > > > ___ > 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 ]
Re: python 3
Hi Nils, FiPy is not written for Python 3, but it does work after applying 2to3. However, only the Scipy solvers are then available. See, https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#running-under-python-3 To fully make the switch to Python 3, FiPy needs to be natively updated, which hasn't happened yet. Also, FiPy needs to move away from using Trilinos as PyTrilinos hasn't been updated to Python 3. Unfortunately, I can't give a timeline on when those changes will be made. Cheers, Daniel On Wed, Jun 7, 2017 at 6:16 AM, Nils Becker <nils.bec...@bioquant.uni-heidelberg.de> wrote: > hi, > > i'm considering updating my scripts from 2.7 to python 3.4 or 3.6. > should i expect any problems? is fipy tested on py 3? > > thanks! > > ___ > 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 ]
Re: solvers for elliptic equations with dominant lower order terms
Hi James, I don't have any hints for preconditioners, but using Trilinos's GMRES and then systematically changing the GMRES parameters and preconditioners might help you find one that works. I have done this in the past, for example, https://gist.github.com/wd15/c28ab796cb3d9781482b01fb67a7ec2d#file-unsegregatedequation-py-L187 You can also get information about how well the solver is converging. Cheers, Daniel On Wed, May 31, 2017 at 3:31 PM, James Pringle <jprin...@unh.edu> wrote: > > Dear all -- > > I need to solve a second order PDE in 2D of the form > > A(x,y)*(eta_xx+eta_yy)+B(x,y)*eta_x+C(x,y)*eta_y = 0 > > This is an advective diffusive balance, with B and C representing the > advective velocities. > > > Over much of the domain the low order terms dominate. For some geometries, > the default solver I am using (SciPy LinearLUSolver) is having trouble find > an answer (especially when the characteristics defined by B and C form > loops. > > Does anyone have recommendations for how to choose a better solver? (I have > blindly tried all the ones in SciPy, the ones besides LinearLUSolver do > worse). > > And are there any hints on how to use preconditioners? There is discussion > of preconditioners and Trilinnos, but I can't find hints on how one would > choose to use them and when they might help. > > Pointers to the broader literature are welcome! > > Thank you, > Jamie Pringle > > ___ > 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 ]
Re: Complex conjugates in FiPY
On Tue, May 23, 2017 at 8:05 AM, Sergio Manzetti <sergio.manze...@fjordforsk.no> wrote: > > Dear Daniel. In this example I was trying to find out a solution for the > complex differential equation: > > i dphi/dt + d²phi/dx^2 + |phi^2|phi = 0 > > However, I just learned from the fipy ring that complex conjugates are not > possible to include in Fipy. Right, FiPy knows nothing about complex numbers. > However, if I should calculate the given eqn without the complex conjugate: > > dphi/dt + d²phi/dx^2 + |phi^2|phi = 0 That is a strange equation because it induces counter diffusion, which probably won't solve. I'm not sure whether the source term does something to counter that, but you can try. > which example should I follow in the manuscript? There isn't an example that matches this equation exactly, but try http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html for the TransientTerm and DiffustionTerm and http://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.simple.html for the ImplicitSourceTerm. > I am not sure what the script does, when one sets a phi value before the > given PDE...when I thought that the phi value was found exactly by FipY? Setting the value before solving the PDEs provide the initial value if it's an initial value problem. Even a non-initial value problem needs to be start somewhere for the non-linear numerical iteration. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Complex conjugates in FiPY
On Sat, May 20, 2017 at 5:02 AM, Sergio Manzetti <sergio.manze...@fjordforsk.no> wrote: > > Dear Daniel, I am wondering if you can clarify a small. thing. > > In the given script, phi is set as e^ix, and the numerical simulation treats > the given PDE. Is phi tested for wether it is a result of the given PDE in > this script ? Or does the script do something else? I'm not quite sure what you're asking, but I'm sure that the script below does not work as you intend it to work. > #!/usr/bin/env python > # testing a non-complex variant of the NLSE > > import numpy > import cmath as math > from fipy import * > from fipy import numerix > > nx = 50 > dx = 1. / float(nx) > > mesh = Grid1D(nx=nx,dx=dx) > X = mesh.cellCenters[0] > > phi = CellVariable(mesh=mesh, name="Solution") > phi.setValue(0.5-0.5*numerix.exp((1j*X))) At this point your script is broken, "phi.value.imag" is all zero while "(0.5-0.5*numerix.exp((1j*X))).value.imag" is non-zero. The type of the CellVariable is wrong initially and the type doesn't change when the value is reset. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Compressible euler equations
On Thu, May 18, 2017 at 2:16 PM, Thibault Bridel-Bertomeu <thibault.bridellel...@gmail.com> wrote: > > 1/ for the pressure equation, when you meant implicitTerm were you thinking > about this : > > eqnP = fipy.ImplicitSourceTerm(coeff=1.0, var=p) == gm1*roE - > 0.5*gm1*(roU**2 + roV**2)/ro Exactly, you can also use implicit terms for roU and roV as well, but linearized. > 2/ Okay, for the dp/dx and dp/dy, I agree now I see that we can use a > centralconvectionterm indeed. Can you confirm they are to be written as : > > coeffForX = fipy.CellVariable(mesh=mesh, rank=1) > coeffForX[0] = 1.0 > coeffForX[1] = 0.0 > fipy.CentralDifferenceConvectionTerm(coeff=coeffForX, var=p) That looks right. > for dp/dx, and as : > > coeffForY = fipy.CellVariable(mesh=mesh, rank=1) > > coeffForY[0] = 0.0 > > coeffForY[1] = 1.0 > > fipy.CentralDifferenceConvectionTerm(coeff=coeffForY, var=p) Seems OK. > > > for dp/dy ? > As for the roE equation, I combined the gradient components into this : > > coeffConv = fipy.CellVariable(mesh=mesh, rank=1) > coeffConv[0] = roU/ro > coeffConv[1] = roV/ro > fipy.CentralDifferenceConvectionTerm(coeff=coeffConv.faceValue, var=p) > > Does it seem reasonable to you ? That last one isn't right. I think that the roU and roV won't be updated as you advance the time steps. I think that you need to define the "coeffConv" to be coeffConv = (fipy.Varialbe((1, 0)) * roU + fipy.Variable((0, 1)) * roV) / ro The assignment into convCoeff that you are using above just inputs the current values from roU / ro into coeffConv not the dependency. > 3/ By non-linear, I get that you mean sweeping the equations several times > before actually updating the solution in time. But I don’t know what > equations i am supposed to be sweeping here, I think I lack the experience. > The whole system should be swept, like so : You need to sweep all of them as you're doing. I don't think that the order matters. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Compressible euler equations
Hi Thibault, The paper and your attempts at solving it are very impressive, nice work. I'm not sure I understand the issue with the equation of state / 2D issue though. Is it that you want to couple in the equation of state with the other equations instead of solving explicitly? Here is a link for some work I did with FiPy quite a long time ago, https://gist.github.com/wd15/c28ab796cb3d9781482b01fb67a7ec2d It resulted in a publication, see https://arxiv.org/pdf/1006.4881.pdf. In that work, everything is being solved in one matrix, see https://gist.github.com/wd15/c28ab796cb3d9781482b01fb67a7ec2d#file-input-py-L245 Look at, https://gist.github.com/wd15/c28ab796cb3d9781482b01fb67a7ec2d#file-input-py-L286 That particular equation is like an equation of state (it has no transient term). The point is that with the use of ImplicitSourceTerm's, you can include an equation of state implicitly. I'm not sure whether this is your issue or not, but it might be related. The above set of equations are dissipative, which yours are not, but also, I wasn't concerned with resolving the shocks numerically so a Roe solver wasn't used. Also, this work was done before FiPy had coupled solutions. If you do include the equation of state in the matrix, then think carefully about how you linearize the "roU" and "roV" terms. Also, as the equation of state is like an immediate update, it might be necessary to relax the updates a bit during the non-linear iteration cycle. Cheers, Daniel On Mon, May 15, 2017 at 11:46 AM, Thibault Bridel-Bertomeu <thibault.bridellel...@gmail.com> wrote: > Hello Daniel, > > Thanks for the answer > > An implementation of the Roe scheme would certainly go a long way in > simplifying the resolution of hyperbolic equations but I think FiPy already > has was it takes to do that - terms wise. > I know of CLAWPACK yes, but I really would like to see if FiPy can handle > fluid dynamics equations. I work at CERFACS, a lab in France centered around > the use of Computational Fluid Dynamics for academic and industrial purpose, > and although we do have a production code, I would like to let my colleagues > know about FiPy because I think it has the potential to be a great > first-approach solver for simple to mildly complex problems. > > I actually already succeeded in solving (and validating with the known Sod > shock tube test case) the 1D compressible Euler equations (script attached > with an exact solver for reference). > But then in 1D, things are pretty easy - equation wise. > If you have time to check the code, can you tell me if you see any kind of > mis-use of FiPy ? Are there any lines you would have written differently ? > > In 2D everything becomes harder and I can’t get it right - although I > suspect it is because I don’t know FiPy very well yet. > If you don’t mind, I have a couple questions to try and make it work. > > 1/ I do not get how to include the resolution of an equation of state in the > system of differential equations. I wrote a stub of code (attached, called > 7-covo2D.py) in which I solve for (rho, rhoU, rhoV, rhoE) as usual. An easy > way to write the Euler equations however is to include a notion of pressure, > which is related to the unknown by an equation of state (it is not a > differential equation though) - see for instance > http://bender.astro.sunysb.edu/hydro_by_example/compressible/Euler.pdf. > Because of this, I have to inject the EoS directly into the equations and > try and make do with that, but the Finite Volume method is not designed to > have the fluxes computed like I do in the script, so naturally everything > crashes or gives nonsensical results. > Do you see a possibility in FiPy to solve 4 differential equations and 1 > algebraic equation (the EoS) in a coupled manner ? > > 2/ If it is not possible to include the equation of state in the system, how > would you take the -dp/dx and the -dp/dy from the momentum equations (with p > expanded as a function of rhoE, rhoU and rhoV) into account ? > > Thank you very much for your help, > > Best regards, > > Thibault > > > > 2017-05-15 15:38 GMT+02:00 Daniel Wheeler <daniel.wheel...@gmail.com>: >> >> Hi Thibault, >> >> I started down the road of implementing a Riemann solver in FiPy, see >> >> >> https://github.com/usnistgov/fipy/blob/riemann/fipy/terms/baseRoeConvectionTerm.py >> >> Unfortunately, I never completed and merged this work. You might want >> to look into CLAWPACK for a code to better handle compressible flow. >> >> I've given a few answers about this in the past as well, see >> >>- http://git.net/ml/python-fipy/2012-02/msg00024.html >> >>- >> http://fipy.nist.narkive.com/A0gJrSl2/semilinear-
Re: Fresh installation tests fail
Hi Carsten,
I'm not sure why that error is occurring, but the following Dockerfile
works for me.
https://github.com/wd15/extremefill2D-dockerize/blob/master/Dockerfile
It might help you debug. The FiPy version is "3.1.3.dev138+gecbe868".
Maybe try that and see if it helps.
Cheers,
Daniel
On Tue, May 16, 2017 at 5:14 PM, Carsten Langrock <langr...@stanford.edu> wrote:
> Following installation guide using conda under macOS 10.12.5, I am not able
> to run the fipy.test() function without it throwing an exception. The output
> is listed below. Any hint on how to solve this would be appreciated.
>
> Python 2.7.13 | packaged by conda-forge | (default, May 2 2017, 13:29:36)
> [GCC 4.2.1 Compatible Apple LLVM 6.1.0 (clang-602.0.53)] on darwin
> Type "help", "copyright", "credits" or "license" for more information.
>>>> import fipy
>>>> fipy.test()
> running egg_info
> creating
> /var/folders/68/q7ztlcxx3nldk6drfnp4l5_hgn/T/tmp6x5SYr/FiPy.egg-info
> writing
> /var/folders/68/q7ztlcxx3nldk6drfnp4l5_hgn/T/tmp6x5SYr/FiPy.egg-info/PKG-INFO
> writing top-level names to
> /var/folders/68/q7ztlcxx3nldk6drfnp4l5_hgn/T/tmp6x5SYr/FiPy.egg-info/top_level.txt
> writing dependency_links to
> /var/folders/68/q7ztlcxx3nldk6drfnp4l5_hgn/T/tmp6x5SYr/FiPy.egg-info/dependency_links.txt
> writing manifest file
> '/var/folders/68/q7ztlcxx3nldk6drfnp4l5_hgn/T/tmp6x5SYr/FiPy.egg-info/SOURCES.txt'
> reading manifest file
> '/var/folders/68/q7ztlcxx3nldk6drfnp4l5_hgn/T/tmp6x5SYr/FiPy.egg-info/SOURCES.txt'
> writing manifest file
> '/var/folders/68/q7ztlcxx3nldk6drfnp4l5_hgn/T/tmp6x5SYr/FiPy.egg-info/SOURCES.txt'
> running test
> running build_ext
> fipy version 3.1.3
> numpy version 1.10.4
> pysparse version 1.3-dev
> scipy version 0.17.1
> matplotlib version 2.0.0
> gist is not installed
> mpi4py version 2.0.0
> Trilinos version: 11.10.2
> PyTrilinos version: 4.10d
> mayavi version 4.5.0
> gmsh version 2.15.0
> Traceback (most recent call last):
> File "", line 1, in
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/site-packages/fipy/__init__.py",
> line 164, in test
> cmdclass={'test': _TestClass(_test)})
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/distutils/core.py",
> line 151, in setup
> dist.run_commands()
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/distutils/dist.py",
> line 953, in run_commands
> self.run_command(cmd)
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/distutils/dist.py",
> line 972, in run_command
> cmd_obj.run()
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/site-packages/setuptools/command/test.py",
> line 211, in run
> self.run_tests()
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/site-packages/fipy/tests/testClass.py",
> line 236, in run_tests
> testLoader = loader_class()
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/unittest/main.py",
> line 94, in __init__
> self.parseArgs(argv)
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/unittest/main.py",
> line 149, in parseArgs
> self.createTests()
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/unittest/main.py",
> line 158, in createTests
> self.module)
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/unittest/loader.py",
> line 130, in loadTestsFromNames
> suites = [self.loadTestsFromName(name, module) for name in names]
> File
> "/Users/langrock/anaconda/envs/MYFIPYENV/lib/python2.7/unittest/loader.py",
> line 100, in loadTestsFromName
> parent, obj = obj, getattr(obj, part)
> AttributeError: 'module' object has no attribute 'testFiPy'
>
> ___
> 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 ]
Re: Compressible euler equations
Hi Thibault, I started down the road of implementing a Riemann solver in FiPy, see https://github.com/usnistgov/fipy/blob/riemann/fipy/terms/baseRoeConvectionTerm.py Unfortunately, I never completed and merged this work. You might want to look into CLAWPACK for a code to better handle compressible flow. I've given a few answers about this in the past as well, see - http://git.net/ml/python-fipy/2012-02/msg00024.html - http://fipy.nist.narkive.com/A0gJrSl2/semilinear-wave-equation-in-fipy - http://fipy.nist.narkive.com/YVZTRM0G/1d-coupled-fluid-equations-in-fipy Cheers, Daniel On Sun, May 14, 2017 at 9:33 AM, Thibault Bridel-Bertomeu <thibault.bridellel...@gmail.com> wrote: > Good afternoon, > > I was wondering if anyone had ever tried implementing the compressible euler > equations in FiPy ? > I have been trying for a while now but even in 1D I can’t get a proper shock > tube solution. > > If anyone has ever tried, I would be infinitely grateful if they were to > share their knowledge !! > > Thank you very much, > > T. BRIDEL-BERTOMEU > > > > ___ > 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 ]
Re: Complex Valued CellVariables?
Hi Mike, There isn't any way to have complex valued fields. Would it be possible to split into real and imaginary parts? Is there a numerical or mathematical reason why that isn't valid? Cheers, Daniel On Tue, May 9, 2017 at 6:02 PM, Michael Waters <waters.mik...@gmail.com> wrote: > Hi, > > I am experimenting with wave equations in FiPy and am wondering if there > is a good way to go about implementing complex valued fields? > > Best, > > -Mike > > ___ > 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 ]
Re: Pulse Laser Propgation Eqn setup
: > return 0 > else: > gamma = ((w*numerix.sqrt(m*Eg))/(e*E)) > gamma1 = (gamma**2)/(1+gamma**2) > gamma2 = 1.0 - gamma1 > x = > (2/numerix.pi)*(Eg/(hbar*w))*(numerix.sqrt(1+gamma**2)/gamma)*ellipe(1/(1+gamma**2)) > Q = numerix.sqrt((numerix.pi)/(2*ellipk(gamma2)))*quad(lambda > j: integ(j,x,gamma1,gamma2),0,1)[0] > keld = > ((2*w)/(9*numerix.pi))*(((w*m)/(hbar*numerix.sqrt(gamma1)))**(3/2))*Q*numerix.exp(-1*numerix.pi*(int(x+1))*((ellipk(gamma1)-ellipe(gamma1))/ellipe(gamma2))) > #keld = keld*100 > return numerix.absolute(keld) > > def multi(E): > tmp=numerix.empty((),) > for i in range(len(E)): > numerix.append(tmp,Keld(E()[i])) > return tmp > > if __name__ == "__main__": > nr = 50 > nz = 50 > dr = 1e-6 > dz = dr > L= dr*nr > mesh = mesh = CylindricalGrid2D(nr=nr, > nz=nz,dr=dr,dz=dz,origin=((-L/2.0,),(0,)))#Grid2D(dx=dx, dy=dy, nx=nx, > ny=ny) > > #create variables > E = CellVariable(name="real part electric Field", mesh=mesh, > hasOld=True, value=0.1, ) > Estar = CellVariable(name="diff real part electric Field", mesh=mesh, > hasOld=True, value=0.0 ) > J = CellVariable(name="imag part electric Field", mesh=mesh, > hasOld=True, value=0.1 ) > Jstar = CellVariable(name="diff imag part electric Field", mesh=mesh, > hasOld=True, value=0.0 ) > rho = CellVariable(name="carrier concentration", mesh=mesh, > hasOld=True, value=1.0e12 ) > time = Variable() > > #make constants for eqn system > #physical constants > c = 299792458.0 > > #material Parameters > n = 2.5 > alpha = 0.05 > n2= 3.5e-16 > vg= 1.0e18 > Beta2 = 1.0e18 > sigma = 1. > tau = 1e-16 > Eg= 3.3*1.602e-19 > #Laser Parameters > wavelength = 1030e-9 #wavelength in nm > w = 2.0*numerix.pi*c/(wavelength) > d = 25e-6 #Focusdepth > wf= 1.0e-6 #Beam radius in Focus > tp= 1e-9 #Puls length > Ein = 100e-3 #Pulsenergy > timeStep = 1e-11 > desiredResidual = 1.0e-8 > totalElapsedTime = 1000*timeStep > > #physical relations > k= n*w/c > > #equation parsmeters > gamma = numerix.array(((1./(2*k), 0.), (0., 0.))) > > #Laserparameter relations > zf = numerix.pi*(wf**2)*n/wavelength > w0 = wf*numerix.sqrt((1.0+((d**2)/(zf**2 > f = d+zf**2.0/d > Pin = Ein/(tp*(numerix.pi/2)) > E0 = numerix.sqrt((2.0*Pin)/(numerix.pi*w0)) > > #set upper values > R,Z = mesh.faceCenters > #mask = ((Y > 19.5e-6)) > E.constrain((E0*numerix.cos((k*(R)**2.0)/(2.0*f)))/(numerix.expR)**2)/(w0**2.0))+((time**2.0)/(tp**2.0, > where=mesh.facesTop) > J.constrain((E0*(-1)*numerix.sin((k*(R)**2.0)/(2.0*f)))/(numerix.expR)**2.0)/(w0**2.0))+((time**2)/(tp**2, > where=mesh.facesTop) > > #setup equation system > eqn1 = TransientTerm(var=E)==Estar > eqn2 = TransientTerm(var=J)==Jstar > eqn3 = E.grad[1] + TransientTerm(coeff=vg, var=E) == > DiffusionTerm(coeff=(-gamma), var=J) \ > +TransientTerm(coeff= Beta2/2.0,var=Jstar)- > (k*n2*((E**2.0))+((J**2.0))*J)+(w0*tau*rho*J)\ > -((sigma*rho*E)/2)-((multi(numerix.sqrt((E**2.0)+(J**2.0)))*E)/(2*((E**2)+(J**2.0 > - alpha*E > eqn4 = J.grad[1] + TransientTerm(coeff=vg, var=J) == > DiffusionTerm(coeff=(gamma), var=E) \ > -TransientTerm(coeff= Beta2/2.0,var=Estar)+ > (k*n2*((E**2))+((J**2))*E)-(w0*tau*rho*E)\ > -((sigma*rho*J)/2)-((multi(numerix.sqrt((E**2.0)+(J**2)))*J)/(2.0*((E**2)+(J**2.0 > - alpha*J > eqn5 = TransientTerm(var=rho)== > rho*thornber((E**2.0)+(J**2.0))+multi((E**2.0)+(J**2.0)) > coupledeqn = eqn1 & eqn2 & eqn3 & eqn4 & eqn5 > > #solve equation system > while time < totalElapsedTime: > residual = 1e5 > E.updateOld() > Estar.updateOld() > J.updateOld() > Jstar.updateOld() > rho.updateOld() > #viewer = Viewer(vars=J) > #viewer.plot() > time.setValue(time() + timeStep) > while residual > desiredResidual: > residual = coupledeqn.sweep(dt=timeStep) > > viewer = Viewer(vars=J) > viewer.plot() > > > ___ > 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 ]
Re: axisymmetric solidification problem in FiPy: trying to start with a basic 1D Stefan problem (solidification)
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=n
Re: derivative boundary condition in 2D problem
On Tue, Jan 24, 2017 at 4:33 AM, Francisco Vega Reyes <fv...@unex.es> wrote: > How do I implement a boundary condition \partial T /\partial y = 0 at > facesTop (for instance)? In that case, I don't believe that you need to add any constraints as that is the natural boundary condition for a diffusion term. It might depend on the equation being solved though in FiPy. > It seems to me that doing var.faceGrad.constrain(((0,),(0,)), > where=mesh.facesTop) (the initial example in the manual) I would be like > doing \partial T /\partial x + \partial T /\partial y = 0... This might be nonsensical, but I think that the constraints have no impact for any quantity that is tangential to the direction of the outbound normal only the normal direction. That could well be a weakness with the notation as we currently have it in 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 ]
Re: Question on accessing internal matrices of the system being solved
Hi Ian, Sorry for the slow response. On Thu, Jan 12, 2017 at 12:20 PM, Campbell, Ian <i.campbel...@imperial.ac.uk> wrote: > > 1) Applying numerix.array() to ‘L’, when ’L’ is of type > 'fipy.matrices.scipyMatrix._ScipyMeshMatrix', creates a zero-dimensional > ndarray, with no shape. This isn’t what we expected because L has diagonal > numerical values & ‘---‘ where its sparse “entries” are. > Our goal is to obtain ‘L’ using your suggested method and then to convert it > into the SciPy sparse.csc_matrix format for further processing. The input to > SciPy’s csc_matrix function must be a rank-2 ndarray, but (reasonably > enough!) this fails when we pass csc_matrix a zero-dimensional ‘L’ matrix. See, https://github.com/usnistgov/fipy/blob/develop/fipy/matrices/scipyMatrix.py#L266 I think you need the "matrix" attribute of "fipy.matrices.scipyMatrix._ScipyMeshMatrix" and I think that is the raw Scipy version of the matrix (whatever format that is). You can then call "toarray()" on that is seems. My previous instructions were wrong. So just using "L.numpyArray" should also achieve the same. > 2) We see from the 2009 paper that it’s a three-point stencil used for the > generation of the discretisation matrix in a first order scheme. What > stencil is used for 2nd order schemes? Depends on the term of course, but for a diffusion term on a square grid it is the same as finite difference which would be a 5 point stencil. The convection terms are mostly first order as currently implemented in FiPy. This is a good book for FV method, http://www.springer.com/us/book/9783319168739, which describes some of the schemes. > 3) How do we implement a higher (e.g. 8th & 12th order central-difference) > order schemes in FiPy? That's not easy at all. I don't think it is designed well enough for that. It would require a major rewrite to easily add new convection schemes. Cheers, Daniel -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Question on accessing internal matrices of the system being solved
Hi Krishna, Yes you can do this, see https://github.com/usnistgov/fipy/blob/develop/fipy/terms/term.py#L332 The following should work, equation.cacheMatrix() equation.solve(...) the_matrix = equation.matrix The cacheMatrix is require because otherwise the reference to the matrix is lost. I'm not exactly certain on the matrix format, but you should be able to do np.arrray(the_matrix) to make it into a 2D numpy array. Hope that helps. On Wed, Jan 11, 2017 at 3:39 PM, Gopalakrishnan, Krishnakumar <krishnaku...@imperial.ac.uk> wrote: > Hi, > > I was wondering if there is a way to access the internal matrices of the > system of equations formulated by Fipy , given the coefficient terms of the > general conservation equation. > > Let's say I have a simple diffusion equation (eg. Heat equation). I can > hand-code the differentiation matrix and BC vector for a constant dx using > the 2nd order central difference formula. > > Is there a way to see the matrix/vector system of the discretised system > that FiPy has internally formed before we call solve() or sweep() ? I am > still learning basics of discretisarion schemes and it will be helpful to > visualise how FiPy internally forms the matrices, and how they compare to > the hand-derived matrices/vectors > > Regards, > > Krishna > > > ___ > 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 ]
Re: Mailing list archive
Hi Ray, Gmain is dead unfortunately. Try this instead, https://www.mail-archive.com/fipy@nist.gov/. Cheers, Daniel On Wed, Nov 23, 2016 at 9:54 AM, Raymond Smith <smit...@mit.edu> wrote: > Hi, FiPy. > > I'm having trouble accessing the mailing list archive, > http://dir.gmane.org/gmane.comp.python.fipy > I can get to that page but I can't search or browse through the archive. I'm > wondering if it's just me or if it's down. > > Ray > > ___ > 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 ]
Re: Adaptive mesh refinement in FiPy?
Hi Anders, No, sorry, there is no adaptive mesh refinement in FiPy. Cheers, Daniel On Mon, Nov 21, 2016 at 2:29 AM, Anders Ericsson <anders.erics...@solid.lth.se> wrote: > Hello, > > > I am wondering if there is any adaptive mesh refinement algorithm > implemented in FiPy? > > If so, is there some sort of manual on how to use it? > > > Thanks and best regards, > > Anders > > > ___ > 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 ]
Re: FW: how to set up data transfer between two adjacent nonuniform meshs
On Tue, Oct 25, 2016 at 12:59 PM, Zhekai Deng <zhekaideng2...@u.northwestern.edu> wrote: > Thanks for the tip. Right now, my code looks like following, where > Outlet_BC_value_new is being updated every time, so does the > Outlet_BC_value. It is my understanding that every timestep, > Outlet_BC_value_new is being over-write with new interpolated value, and > being assign to the Variable, which automatically update it with the phi > constrain. Is this somewhat similar to what you are suggesting for the cache > the memory space ? Yes, however, profile first to see if it is relevant in your case. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FW: how to set up data transfer between two adjacent nonuniform meshs
On Tue, Oct 25, 2016 at 12:34 PM, Zhekai Deng <zhekaideng2...@u.northwestern.edu> wrote: > Thank you. I implemented as you suggested and it works very good now. This > approach seems also applicable for the mesh that are not aligned in the > interface. Exactly. One thing though, it is probably a good idea to cache the mapping from one grid to the other as this is expensive to calculate (especially of you do this at every time step or sweep). We provided this functionality in FiPY, I'm not sure if Scipy lets you do that. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FW: how to set up data transfer between two adjacent nonuniform meshs
On Mon, Oct 24, 2016 at 12:20 AM, Zhekai Deng <zhekaideng2...@u.northwestern.edu> wrote: > > This is not a issue if mesh system is 1D, and boundary cell is just 1. Or > maybe on the rectangular mesh system. However, does any one has any idea to > how to transfer it properly ? Hi Zhekai; I haven't read through everything but I would guess that you need to interpolate from one mesh to another. To do that you can use one of Scipy's inerpolation functions such as https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html We also built that capability into FiPy, but it's probably better to use Scipy's interpolation at this point. https://github.com/usnistgov/fipy/blob/develop/fipy/variables/cellVariable.py#L167 If the mesh points are aligned, it is still a good idea to use the interpolation approach as it removes the issue of worrying about mesh ordering. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Solving convection only problem
Hi Zhekai, There is generally a lot of numerical diffusion when solving convection problems with first order schemes and even some numerical diffusion when using higher order schemes. There are many different schemes and a mass of literature on how to preserve square waves, shocks and hyperbolic equations, but FiPy doesn't have any of those schemes implemented (e.g. TVD schemes come to mind). There is also a secondary issue when coupling hyperbolic equations to do with how the flux is calculated that FiPy doesn't address (the Riemann problem, roe solvers etc). However, for many convection-diffusion problems the time scale of the shocks is not worth resolving or is impossible to resolve while also resolving much longer time scales. When resolving at the convection time scale there is often no benefit from the implicit schemes that FiPy uses. Basically, FiPy is not a great tool for shock problems. CLAWPACK may be something that you could look at for this. I think it's fully explicit and it's focus is on hyperbolic coupled equations. I hope that helps. Cheers, Daniel On Thu, Oct 20, 2016 at 12:58 AM, Zhekai Deng <zhekaideng2...@u.northwestern.edu> wrote: > Hi all, > > I am trying to use Fipy to solve convection only problem for the > concentration moved only by solid body rotation in a "circular" shape > geometry. > > By looking at the examples online, I found out that > > http://www.ctcms.nist.gov/fipy/examples/convection/generated/examples.convection.source.html#module-examples.convection.source > > and some of the level set example appears to allow me do it. > > I implemented the approach from convection example. However, the solution > still looks has diffusion ( or maybe artificial smoothness ) as the > concentration move with velocity field. I have attached my example code, > that concentration enter from the top right side of the geometry, and > undergo solid body rotation eventually to the left side, and flow out of the > domain. So my question is that is there any way to further reduce the > diffusion? Also, does anyone know where this "diffusion" is coming from ? > > The approach that I have tried but did not work are following: > 1. Solve the equation with very small diffusion coefficient (1e-8) > 2. Reduce the timestep or refine the mesh size does not seem to help very > much > > I have attached my example code in this email. Thank you very much. > > Best, > > Zhekai > > > ___ > 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 ]
Re: Question on time-marching scheme in Fipy
Sorry for the slow response, see answers below. On Mon, Oct 10, 2016 at 6:08 PM, Gopalakrishnan, Krishnakumar <krishnaku...@imperial.ac.uk> wrote: > > However, my questions are more general, to be executed when updateOld() is > called. > > · What’s the default implicit scheme in fipy? It's fully implicit unless the user sets it up in a different way. > · How does one go about implementing a specific 2nd order > time-stepping scheme such as (Adams-Bashforth, BDF etc.) We don't have any easy way to do higher order time stepping right now and I'm not aware of any attempts to do so. > · Is there any way to use the FVM only for the spatial > discretisation, i.e. use a method of lines approach for the time-stepping ? It might be possible with source terms, but I haven't tried. > I apologise if the questions sound too basic here. I am just curious about > understanding fipy’s default scheme and implementing an own time-stepper. It's a very good question, but I don't have any helpful answers. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Memory Leakage & Object Build-up with FiPy Sweeps
Hi Ian, Could you possible post your code or a version of the code that demonstrates the problem? Also, do you have the same issue with different solver suites? Cheers, Daniel On Fri, Sep 30, 2016 at 12:41 PM, Campbell, Ian <i.campbel...@imperial.ac.uk > wrote: > Hi All, > > > > We are sweeping six PDEs in a time-stepping loop. We’ve noticed that as > CPU time progresses, the duration of each time-step increases, although the > sweep count remains constant. This is illustrated in the Excel file of data > logged from the simulation, which is available at the first hyperlink below. > > > > Hence, we suspected a memory leak may be occurring. After conducting > memory-focused line-profiling with the vprof tool, we observed a linear > increase in total memory consumption at a rate of approximately 3 MB per > timestep loop. This is evident in the graph at the second link below, which > illustrates the memory increase over three seconds of simulation. > > > > As a further step, we used Pympler to investigate the source of RAM > consumption increase for each timestep. The table below is an output from > Pympler’s SummaryTracker().print_diff(), which describe the additional > objects created within every time-step. Clearly, there are ~3.2 MB of > additional data being generated with every loop – this correlates perfectly > with the total rate of increase of memory consumption reported by vprof. > Although we are not yet sure, we suspect that the increasing time spent per > loop is the result of this apparent memory leak. > > > > We suspect this is the result of the calls to .sweep, since we are not > explicitly creating these objects. Can the origin of these objects be > traced, and furthermore, is there a way to avoid re-creating them and > consuming more memory with every loop? Without some method of unloading or > preventing this object build-up, it isn’t feasible to run our simulation > for long durations. > > dict > > 2684 > > 927.95 > > KB > > type > > 1716 > > 757.45 > > KB > > tuple > > 9504 > > 351.31 > > KB > > list > > 4781 > > 227.09 > > KB > > str > > 2582 > > 210.7 > > KB > > numpy.ndarray > > 396 > > 146.78 > > KB > > cell > > 3916 > > 107.08 > > KB > > property > > 2288 > > 98.31 > > KB > > weakref > > 2287 > > 98.27 > > KB > > function (getName) > > 1144 > > 67.03 > > KB > > function (getRank) > > 1144 > > 67.03 > > KB > > function (_calcValue_) > > 1144 > > 67.03 > > KB > > function (__init__) > > 1144 > > 67.03 > > KB > > function (_getRepresentation) > > 1012 > > 59.3 > > KB > > function (__setitem__) > > 572 > > 33.52 > > KB > > SUM > > 3285.88 > > KB > > > > > > https://imperialcollegelondon.box.com/s/zp9jj67du3mxdcfgbc4el8cqpxwnv0y4 > > > > https://imperialcollegelondon.box.com/s/ict9tnswqk9z57ovx8r3ll5po5ccrib9 > > > > With best regards, > > > > - Ian & Krishna > > > > P.S. Daniel, thank you very much for the excellent example solution you > provided in response to our question on obtaining the sharp discontinuity. > > > > Ian Campbell | PhD Candidate > > Electrochemical Science & Engineering Group > > Imperial College London, SW7 2AZ, United Kingdom > > > > ___ > 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 ]
Re: FiPy sweep convergence bottoms out
> >> >> >> How do we get our residuals to drop faster, i.e. with lesser sweeps and >> faster convergence ? I tried changing solvers and tolerances, but curiously >> enough the results remain identical. >> >> >> >> Any pointers on this will be much appreciated. >> >> >> >> >> >> Krishna >> >> >> >> >> >> >> >> >> ___ >> 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 ] > > > > > ___ > 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 ]
Re: sources errors in advection/diffusion problems, and one solution
is initially set to. If the line > >psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0) > >is changed so psi is initially 1, the error goes away entirely; if >it is set to some other value, you get different errors. I do not >know if this is a bug, or just numerical error in a well programmed >solver. > >Now if you run SquareGrid_HomemadeDelaunay_transient which implements > > psi_t = \Del\dot (D\del psi + u*psi) > >you can see that the error in the numerical solution is advected >out of the x=0 boundary, and the solution approaches the true >solution of psi=1 rapidly. > >The interesting question is if the formulation of the boundary >condition at x=0 could be altered to less excite the spurious mode? > >Also, why does the "initial" condition have any effect on the >steady equation? > >Cheers, >Jamie > > > ___ > 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 ]
Re: JOR solver to help in accelerating sweeps?
Hi Krishna, Have you tried thoroughly profiling your code? You may want to check how much of the time is being spent making copies of arrays and array arithmetic as opposed to the linear solve before changing the solvers for efficiency purposes. I'm not sure the linear solver matters for reducing the number of sweeps. It may require a better non-linear solver and using a single matrix. In the long run, it might be worth contemplating some of the non-linear solvers in Trilinos that remove the boundary between linear and non-linear solvers and also use matrix free methods. Sorry that I can't be more helpful. On Thu, Sep 15, 2016 at 11:06 AM, Gopalakrishnan, Krishnakumar <krishnaku...@imperial.ac.uk> wrote: > > Hello, > > > > We have a specific question about how a solver helps in accelerating > convergence of a loosely coupled system of PDEs. > > > > Due to the fundamentally different properties/behaviour/geometry, we have a > loosely coupled system of 8 PDEs defined along 5 different meshes, with data > manually copied from certain regions of one mesh to suitable regions of > another after sweeping each PDE. > > > > Quite unsurprisingly, the system of equations needs a lot of sweeps (80 to > 100) within each time-step to converge, and consequently the simulation model > isn’t useable. Although, the trilinos based parallel solvers are correctly > set-up and working correctly for the example problems shipped with FiPy, the > mpirun command complains when invoked on our (complicated script). We assume > that it is due to the complicated setup of meshes and the classes/methods and > objects that we have instantiated for our application. The serial code runs > just fine. > > > > So, in order to speed up our simulation, we looked to implement dynamic > relaxation factors for the sweeps. This has only a limited success so far. > > > > Next, I wonder whether including a Successive over-relaxation solver will > help in this case ? There are 8 PDEs to be swept continuously until the > residue of all of them gets lowered below a tolerance. Does the Fipy’s JOR > solver (wrapper for Pysparse Jacobi/over-relaxation solvers) help in this > case ? is there any other way to speed up our system ? > > > > > > Best Regards > > > > Krishna > > > > > > > > > > > _______ > 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 ]
Re: Interesting (perhaps unique) approach to parallelising FiPy sweeps
On Fri, Sep 16, 2016 at 3:38 PM, Gopalakrishnan, Krishnakumar <krishnaku...@imperial.ac.uk> wrote: > > The questions are : 1) Has anybody tried this before ? 2) Does this sound > problematic from a technical or practical perspective. 3) Does this sound > remotely useful. To my knowledge, this hasn't been done with FiPy. Certainly, this can be set up with FiPy using mpi4py or ipython's parallel tools, which you seem to have done below. It could well be useful in specific cases, but In general this approach is very limited since you can't scale to different number of processors so I would definitely not recommend making too much effort. > I have successfully implemented a multiprocessor parallelisation for > examples.diffusion.coupled (shown below) . The values of the simulation > results (for v0 & v1) after 100 time-steps are very very close to that of the > serial sweep. Although this code is several times slower than the serial > code, we tried this more as a concept demonstrator before embarking to spend > significant effort to convert our serial code to multiprocessor approach. If > you see a red flag, it would be much appreciated if you can help us by > pointing it out. I don't see a red flag there. Try to set it up so that the equations are not redefined for every sweep. Also, do an in place update of the variables using [:] or setValue (return the array rather than the variable from get_res). -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Solving advection, segregation and diffusion equation
On Wed, Sep 14, 2016 at 7:52 PM, Guyer, Jonathan E. Dr. (Fed) <jonathan.gu...@nist.gov> wrote: > > I don't know offhand. With a Péclet number of 100, your problem is almost > completely hyperbolic, which FiPy (and cell-centered Finite Volume) isn't > very good at. Daniel knows more about this and may have some suggestions. > > You might consider adding a TransientTerm for artificial relaxation and > trying the VanLeerConvectionTerm. Definitely use a TransientTerm and see the time step you can get away with. Don't sweep, just use time steps. Using the VanLeerConvetionTerm isn't necessary as the accuracy is not important in a steady state 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 ]
Re: Coupling Boundary Condition of one PDE with source term of another PDE
On Tue, Sep 13, 2016 at 12:59 PM, Guyer, Jonathan E. Dr. (Fed) <jonathan.gu...@nist.gov> wrote: > I'm not sure; Daniel might have a better idea. > > One thing you can do to manage this (if the anisotropic diffusion is > necessary) is change 2c) to: > >c) eq1 = fp.DiffusionTerm(coeff=S * [[[1., 0], [0., 0.]]]) == f(w, a) > # matrix co-oeff has now been simplified by using a single facevariable whose > values are appropriately zeroed out in the 2D grid Krishna, the above is correct the coefficient can be both a rank 2 cell (or face) variable and have it's value vary in space. You could do my_coeff = fp.CellVariable(mesh=mesh, rank=2, value=0.) and then set the value explicitly my_coeff[0, 0] = x**2 ## or however you want to set the values and then set the coeff in the diffusion term to be "my_coeff". This is exactly equivalent to what Jon is saying and similar to what you suggested. It should work fine. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Periodic BC
On Wed, Sep 14, 2016 at 4:12 AM, Francisco Vega Reyes <fv...@unex.es> wrote: > > Hi again, > > is there a way to implement periodic boundry conditions in fipy? Yes, it's done with the mesh class rather than a constraint. A few examples, http://stackoverflow.com/questions/36177369/fipy-simple-convection https://github.com/usnistgov/fipy/blob/64f7866c8dbe50e2c36adfd23c464a53e4a2b763/examples/diffusion/steadyState/mesh1D/inputPeriodic.py Hope that helps. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Extremely Slow .sweep & Updating CellVariable-dependent Boundary Conditions
Hi Ian, Sorry for the slow response. On Wed, Sep 7, 2016 at 5:57 PM, Campbell, Ian <i.campbel...@imperial.ac.uk> wrote: > > Firstly, the FiPy code I’m working with sweeps for four inter-dependent > CellVariables / PDEs in a loop until their residuals falls below specified > tolerances. After profiling the code, it’s clear that this sweeping process > accounts for ~95% of the time spent solving in a given timestep. > Consequently, advancing in time is painfully slow. Assuming that the > CellVariable tolerances are fixed at their largest allowable values, what > would be the best strategy for me to take to reduce the time spent sweeping > these variables? I'm surprised that it's only 95%, it should be closer to 100% as that's mostly all that happens in any given time step. The sweep consists of the matrix build and the linear solve. That's the proportion that you want to understand. Best to profile to see how long the linear solution is as a proportion of the overall sweep time and work from there. Then you can then focus your attention on improving either the linear solve or improving the coefficient calculations (using --inline for instance) depending on what's slow or see if there is a calculation (e.g. whether you use numpy.maximum or max can make a huge difference) that is spending the bulk of the time. > Secondly, some (but not all) of the boundary conditions for these four PDEs > need to be updated between sweeping one CellVariable and the next. Notably, > this is within a single timestep because the boundary conditions are > functions of the CellVariables being swept. What is the recommended approach > for updating the values of these boundary conditions, and why? If the boundary conditions values are defined as variables then they should automatically update themselves or you can update them during the sweep manually with some code if you don't trust the lazy evaluation. Either way that is an explicit approach to updates and will probably have some limitation in the time step. If you can build implicitness into the boundary conditions then that could save you on the time step restriction, but this is highly problem dependent. I think it's best to have a concrete example. > At present, > the following is the approach used: define the boundary condition only > within the loop using x.constrain(value, where=some_boundary). i.e. the > boundary condition is re-defined with new variables at every pass through > the loop. My feeling is that if "value" is an expression involving other variables then it will update automatically so you won't have to keep redefining the constraint during the loop. I think whatever approach you use it shouldn't be necessary to redefine the same constraint over and over. There could be a lot wrong with that approach. You can always make "value" a variable and update it's value explicitly in the loop rather than redefining the constraint. > However, I suspect that this approach may be incorrect, and that a better > approach may be to define the boundary conditions only once, outside of the > loop, and to update them by employing the .value approach for CellVariables. > I am at a loss to understand why this approach may be better, though. Yup. It's just that you might be defining many constraints for a CellVariable. Are they updated or do they all apply or what? I'm not sure, but I don't think they were designed to work like that. Have you tried defining the constraint once and updating the variable's value? Does it give the same results? Cheers, Daniel -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Dynamic under-relaxation factors for FiPy sweep
On Mon, Aug 29, 2016 at 2:40 PM, Gopalakrishnan, Krishnakumar <krishnaku...@imperial.ac.uk> wrote: > > I am trying to sweep for 8 field variables in FiPy, until all their residues > die down to appreciably small values. > > > > Currently, the convergence of this iterative loop is very slow. I have fixed > under-relaxation implemented right now, but looking for implementing a > dynamic under-relaxation factor for each of these 8 variables – something > like an Aitken’s or least-squares method should be helpful. > > > > Is there any helpful resource that the users of FiPy may be able to point me > to ? I don't think there are any FiPy specific resources on that. Your issue is with the algorithm or theory rather than how to implement it in Python / FiPy, right? -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: reaction-advection-diffusion; problems making the solution converge
Hi Nils,
It seems like the code is converging well as far as I can tell. For
example, when I modify the sweep loop to be
print
print 'step: {0}'.format(tt)
while res > 1e-9: #Solve.sweeps:
res = eq.sweep(var=rho, dt=Solve.step,
solver=Solve.customSolver)
print res
every step seems to converge in 5 steps to 1e-10.
step: 20.0
0.00192315625797
4.5626520206e-05
1.08935792171e-06
2.60975675747e-08
6.26480480874e-10
step: 21.0
0.00195977695607
4.651749367e-05
1.11103436368e-06
2.66255569394e-08
6.3934852818e-10
step: 22.0
0.00199768415442
4.74384979099e-05
1.13342888153e-06
2.71708361386e-08
6.52634229319e-10
The non-linear residual only really tells you how well each step
individually is converging. It seems like you were only using the very
first non-linear residual, which uses the old value for the
calculation, which doesn't tell you anything useful. Basically "A_new
* x_old -b_new" is not a measure of anything useful. It is only useful
as a quantity for normalization. FiPy uses the old value or the
previous sweep value to calculate the non-linear residual. So the
residual is
A(n) * x(n-1) - b(n)
where n is the sweep number, and x_(-1) is x_old if the first sweep is
n=0. If only one sweep is executed then the residual isn't a useful
quantity. The linear residual is not returned from the "sweep" method,
which is "A(n) * x(n) - b(n)" and is assumed to be as accurate as the
tolerance of the linear solve. Basically, collecting the first
non-linear residual at each time step is not a useful exercise to
generate a metric for convergence across time steps. In fact, I'm not
even sure how to do that or what it means. I only know how to
understand convergence for a given time step.
One thing that I'm confused about is that it requires 5 sweeps to
converge for a fully linear equation. Is it actually linear? I
couldn't tell from the code.
Cheers,
Daniel
On Fri, Aug 5, 2016 at 8:30 AM, Nils Becker
<nils.bec...@bioquant.uni-heidelberg.de> wrote:
>>> I'm not sure if using ||lhs|| is a good way to normalize the residual.
>>> Surely that could be very small. Is ||lhs - rhs|| getting smaller
>>> without normalization?
>
> if you don't see any obvious errors in building the solution in fipy, i
> could continue poking around if i knew how to evaluate the various terms
> of the solution and their errors at different times. i see methods
> eq_diffprolif.justErrorvector, but i don't quite understand what exactly
> this returns. is it lhs-rhs? what is justResidualVector? can i get those
> for the individual terms on the rhs as well?
>
> cheers, n.
>
> ___
> 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 ]
Re: reaction-advection-diffusion; problems making the solution converge
On Wed, Aug 3, 2016 at 9:27 AM, Nils Becker <nils.bec...@bioquant.uni-heidelberg.de> wrote: >> Could you possible post your code so I can play around with it? I >> can't think offhand what could be the issue or exactly what the error >> is that you're calculating, it seems like the normalized residual or >> something. > > hi daniel, thanks for your reply. the way i was calculating the residual > is to export the solution on a grid, import that into mathematica, > interpolate in space and time and have mathematica calculate > ||lhs-rhs|| / ||lhs|| where ||...|| is the L2 norm. what i saw was that > the relative error was low in the beginning but increased over time, > especially in places where the solution starts to grow exponentially. I'm not sure if using ||lhs|| is a good way to normalize the residual. Surely that could be very small. Is ||lhs - rhs|| getting smaller without normalization? -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Applying Bounding Limits to PDE Solution Variables
On Tue, Aug 2, 2016 at 2:02 PM, Campbell, Ian <i.campbel...@imperial.ac.uk> wrote: > > I've tried the mean-square and RMS values, but they're both orders of > magnitude away from the residual returned by '.sweep'. As such, can you > please tell me the formula used by '.sweep' to compute the residual, so that > I may use the same one? I think it's the L2 norm of "A * x - b" without any normalization, where x is the variable value before the linear solve. See, https://github.com/usnistgov/fipy/blob/develop/fipy/solvers/solver.py#L135 -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: reaction-advection-diffusion; problems making the solution converge
Hi Nils, Could you possible post your code so I can play around with it? I can't think offhand what could be the issue or exactly what the error is that you're calculating, it seems like the normalized residual or something. This is a linear equation so I'm not entirely sure why the non-linear residual wouldn't be very close to zero if the linear solver is converging. Probably best just to see the code rather than speculate. Cheers, Daniel On Mon, Aug 1, 2016 at 7:50 AM, Nils Becker <nils.bec...@bioquant.uni-heidelberg.de> wrote: > hi, > > i scanned the list and the documentation but have not been able to make > this work properly yet: > > i have a 2D pde of the form > > \[ > > \partial_t \rho = \nabla \cdot (\rho \nabla \phi) + D \nabla^2 \rho + > \lambda \rho > > \] > > where phi is a potential field and lambda is a source field, both are > constant in time and smoothly varying in space. > > i have been solving this with fipy, using > > DiffusionTermCorrection(coeff=[D,]) for the diffusion > PowerLawConvectionTerm(coeff=v ) where v = -phi.grad for the advection > ImplicitSourceTerm(coeff=lambda) for the source > > the boundary conditions are either no-flux or dirichlet. i think > advection is dominant, although i can't say i have actually tried to > quantify that. > > i am testing the solution by subtracting the rhs and lhs of the equation > and comparing that to \partial_t \rho in magnitude. i am having trouble > making this error smaller than around 10%: at points where the density > grows strongly, the solution deviates. > > at the moment i am using a loop like this: > > for (st, stnext) in zip(Solve.sample_t, Solve.sample_t[1:]): > for tt in nx.arange(st, stnext, Solve.step): > rho.updateOld() > res = 1. > while res > Solve.tolerance: > res = eq.sweep(var=rho, dt=Solve.step) > > i've played around with making the mesh smaller -- did not seem to make > any difference; reducing the tolerance -- some effect, and reducing the > temporal step size, which did not help much. > > should i be doing something differently? > > thanks for any help! > > n. > _______ > 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 ]
Re: error in gradient even with orthogonal grid
Hi Jamie,
Sorry for the delayed reply.
On Thu, Jul 21, 2016 at 11:31 AM, James Pringle <jprin...@unh.edu> wrote:
>
> However, I am still getting large errors in estimates of the gradient even
> when the mesh is nearly perfectly orthogonal. These errors DO NOT become
> smaller as resolution is increased. I have two question:
>
> Are these errors to be expected in unstructured meshes?
Cell centered FV has two issues with unstructured meshes. You know
about the orthogonality issue, but there is a second issue. The
docstring from the following code might help explain,
https://github.com/usnistgov/fipy/blob/develop/fipy/variables/cellVariable.py#L257
See, `\frac{1}{V_P} \sum_f \vec{n} \phi_f A_f`. This is the
discretized gradient calculation. In this formulation using the
triangular mesh, the calculated \phi_f is not the correct average for
the face as the line between the cell centers doesn't cross the face
in the middle of the face. I think this is called non-conjunctionality
in finite volume speak. See equation 9.8 in the following link for
further explanation (we should probably implement this approach).
http://bit.ly/2aqRiHt
I believe that this is the source of the error in the gradient
calculation. It's a consistent error independent of the mesh
refinement I think (or maybe first order). An alternative gradient
calculation is the least squares gradient, see
https://github.com/usnistgov/fipy/blob/develop/fipy/variables/cellVariable.py#L270
Using that
import fipy as fp
nx = 25
dx = 1. / nx
mesh = fp.Tri2D(nx=nx, ny=nx, dx=dx, dy=dx)
psi = fp.CellVariable(mesh=mesh)
psi.constrain(1.0, where=mesh.facesLeft)
psi.faceGrad.constrain(1.0, where=mesh.facesRight)
eq = fp.DiffusionTerm().solve(psi)
print psi.leastSquaresGrad[:, :100]
seems to give lots of 1s and 0s which seems correct. The cells next to
the boundary are not correct as they aren't using the boundary
conditions in the calculation. I'm not sure about the order of
accuracy of the least squares or why it happens to work for this mesh.
> Is it acceptable to average the gradient over space to reduce the error?
I don't think so unless it can be shown to improve the accuracy with analysis.
--
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Improving non-orthogonality issue in FiPy
Hi Michael, I think this paper explains the issue better than I can, https://www.researchgate.net/publication/263972969_ANALYSIS_OF_THE_NON-ORTHOGONALITY_CORRECTION_OF_FINITE_VOLUME_DISCRETIZATION_ON_UNSTRUCTURED_MESHES Figure 2 shows the control volume and Figure 3 shows the skewness issue. Cheers, Daniel On Wed, Jul 20, 2016 at 4:30 PM, Michael Waters <waters.mik...@gmail.com> wrote: > Hi, I'm not sure I understand the issue. Maybe you could explain more? I > made a diagram to help. > > Cheers, > > -Mike Waters > > > > On 07/20/2016 02:51 PM, Daniel Wheeler wrote: >> >> FiPy has issues with non-orthogonal meshes (i.e. when the vector >> between the cell center isn't parallel to the face normal). We did >> make an attempt at fixing this, which resulted in two diffusion terms, >> >> >> https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermCorrection.py >> >> and >> >> >> https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermNoCorrection.py >> >> As I recall, the "corrected" diffusion term was quite unstable for >> moderate to high non-orthogonality. >> >> I'd like to address this issue again and I was wondering if anyone had >> opinions regarding the latest / best / most tractable approach for >> dealing with non-orthogonality in CC-FV. I haven't followed the >> literature on this for quite some time. >> > > > ___________ > 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 ]
Re: Accuracy of DiffusionTerm for non-uniform mesh
On Wed, Jul 20, 2016 at 5:45 PM, Raymond Smith <smit...@mit.edu> wrote:
> That makes sense. So, here
> https://gist.github.com/raybsmith/b0b6ee7c90efdcc35d6a0658319f1a01
> I've changed it so that the error is calculated as
> sum( (\phi-phi*)^2 * mesh.cellVolumes )
> and depending on the value I choose for the ratio of dx_{i+1} / dx_i, I get
> convergence for exponential spacing ranging from 2nd order (at that ratio =
> 1) and decreasing order as that ratio decreases from unity.
Thanks for that. I'm not entirely sure, but the integral as written on
line 44,
https://gist.github.com/raybsmith/b0b6ee7c90efdcc35d6a0658319f1a01#file-fipy_accuracy-py-L44
may still only be first order accurate for non-uniform grids. I'll try
and investigate this though and see if it matters.
--
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Accuracy of DiffusionTerm for non-uniform mesh
On Wed, Jul 20, 2016 at 1:30 PM, Raymond Smith <smit...@mit.edu> wrote: > Hi, FiPy. > > I was looking over the diffusion term documentation, > http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#diffusion-term > and I was wondering, do we lose second order spatial accuracy as soon as we > introduce any non-uniform spacing (anywhere) into our mesh? I think the > equation right after (3) for the normal component of the flux is only second > order if the face is half-way between cell centers. If this does lead to > loss of second order accuracy, is there a standard way to retain 2nd order > accuracy for non-uniform meshes? This is a different issue than the non-orthogonality issue, my mistake in the previous reply. > I was playing around with this question here: > https://gist.github.com/raybsmith/e57f6f4739e24ff9c97039ad573a3621 > with output attached, and I couldn't explain why I got the trends I saw. > The goal was to look at convergence -- using various meshes -- of a simple > diffusion equation with a solution both analytical and non-trivial, so I > picked a case in which the transport coefficient varies with position such > that the solution variable is an arcsinh(x). I used three different styles > of mesh spacing: > * When I use a uniform mesh, I see second order convergence, as I'd expect. > * When I use a non-uniform mesh with three segments and different dx in each > segment, I still see 2nd order convergence. In my experience, even having a > single mesh point with 1st order accuracy can drop the overall accuracy of > the solution, but I'm not seeing that here. > * When I use a mesh with exponentially decreasing dx (dx_i = 0.96^i * dx0), > I see 0.5-order convergence. That's strange. Are you sure that all the normalization for grid spacing is correct when calculation the norms in that last case? > I can't really explain either of the non-uniform mesh cases, and was curious > if anyone here had some insight. I don't have any immediate insight, but certainly needs to addressed. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Improving non-orthogonality issue in FiPy
FiPy has issues with non-orthogonal meshes (i.e. when the vector between the cell center isn't parallel to the face normal). We did make an attempt at fixing this, which resulted in two diffusion terms, https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermCorrection.py and https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermNoCorrection.py As I recall, the "corrected" diffusion term was quite unstable for moderate to high non-orthogonality. I'd like to address this issue again and I was wondering if anyone had opinions regarding the latest / best / most tractable approach for dealing with non-orthogonality in CC-FV. I haven't followed the literature on this for quite some time. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Accuracy of DiffusionTerm for non-uniform mesh
Hi Raymond,
I was just corresponding with James Pringle off the list about the
same issue for a Delaunay triangulated mesh. Here is my explanation to
him.
I think this is the non-orthogonality issue. I should have realized
this much sooner. The only place we seem to reference this issue is in
www.ctcms.nist.gov/fipy/documentation/numerical/discret.html where it
is stated "this estimate relies on the orthogonality of the mesh, and
becomes increasingly inaccurate as the non-orthogonality increases.
Correction terms have been derived to improve this error but are not
currently included in FiPy [13]."
See this for example,
https://www.researchgate.net/publication/242349760_Finite_volume_method_for_the_solution_of_flow_on_distorted_meshes
In the test mesh that you show, some of the triangles are highly
non-orthogonal to each other. This results in errors. You can view the
quantity of non-orthogonality in the mesh using
nonorth_var = fipy.CellVariable(mesh=mesh, value=mesh._nonOrthogonality)
fipy.Viewer(nonorth_var).plot()
The following code uses triangles, but with a perfectly orthogonality
and, thus, gives perfect results.
import fipy as fp
nx = 50
dx = 1. / nx
mesh = fp.Tri2D(nx=nx, ny=nx, dx=dx, dy=dx)
psi = fp.CellVariable(mesh=mesh)
psi.constrain(1.0, where=mesh.facesLeft)
psi.faceGrad.constrain(1.0, where=mesh.facesRight)
eq = fp.DiffusionTerm().solve(psi)
fp.Viewer(psi - mesh.cellCenters[0]).plot()
nonorth_var = fp.CellVariable(mesh=mesh, value=mesh._nonOrthogonality)
fp.Viewer(nonorth_var).plot()
raw_input('stopped')
As far as your Delaunay triangulation is concerned, finding a way to
smooth out / reduce the non-orthogonality would be one way to improve
the accuracy.
The FiPy docs definitely need to be more explicit about this issue.
Also, actually addressing the orthogonality issue for diffusion terms
at least would be good.
Cheers,
Daniel
On Wed, Jul 20, 2016 at 1:30 PM, Raymond Smith <smit...@mit.edu> wrote:
> Hi, FiPy.
>
> I was looking over the diffusion term documentation,
> http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#diffusion-term
> and I was wondering, do we lose second order spatial accuracy as soon as we
> introduce any non-uniform spacing (anywhere) into our mesh? I think the
> equation right after (3) for the normal component of the flux is only second
> order if the face is half-way between cell centers. If this does lead to
> loss of second order accuracy, is there a standard way to retain 2nd order
> accuracy for non-uniform meshes?
>
> I was playing around with this question here:
> https://gist.github.com/raybsmith/e57f6f4739e24ff9c97039ad573a3621
> with output attached, and I couldn't explain why I got the trends I saw.
> The goal was to look at convergence -- using various meshes -- of a simple
> diffusion equation with a solution both analytical and non-trivial, so I
> picked a case in which the transport coefficient varies with position such
> that the solution variable is an arcsinh(x). I used three different styles
> of mesh spacing:
> * When I use a uniform mesh, I see second order convergence, as I'd expect.
> * When I use a non-uniform mesh with three segments and different dx in each
> segment, I still see 2nd order convergence. In my experience, even having a
> single mesh point with 1st order accuracy can drop the overall accuracy of
> the solution, but I'm not seeing that here.
> * When I use a mesh with exponentially decreasing dx (dx_i = 0.96^i * dx0),
> I see 0.5-order convergence.
>
> I can't really explain either of the non-uniform mesh cases, and was curious
> if anyone here had some insight.
>
> Thanks,
> Ray
>
> ___
> 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 ]
Re: Diffusion-Convection-Reactive Chemisty
Hi Daniel,
It is giving a different error for me
Traceback (most recent call last):
File "ConversionModel.py", line 78, in
V.constrain(V0,mesh.facesLeft)
NameError: name 'V' is not defined
Is this the correct script?
On Wed, Jul 13, 2016 at 3:27 PM, Daniel DeSantis <desan...@gmail.com> wrote:
> Hello,
>
> I am having some trouble getting a workable convection coefficient on a
> reactive chemistry model with a vector velocity. I have reviewed several
> examples from the FiPy mailing list and tried several of the variations
> listed there. Originally, I was receiving the error that a coefficient must
> be a vector. Now I seem to be getting an error that says the index is out of
> bounds. Specific traceback is shown below, along with the code attachment.
>
> Could anyone suggest what I am doing wrong and how to fix it?
>
> I'm relatively new to Python and FiPy, specifically, so please bear with me
> as I am learning.
>
> Thank you all!
>
> --
> Daniel DeSantis
>
> Traceback (most recent call last):
>
> File "", line 1, in
> runfile('C:/Users/ddesantis/Dropbox/PythonScripts/CFD
> Models/ConversionModel.py',
> wdir='C:/Users/ddesantis/Dropbox/PythonScripts/CFD Models')
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py",
> line 699, in runfile
> execfile(filename, namespace)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py",
> line 74, in execfile
> exec(compile(scripttext, filename, 'exec'), glob, loc)
>
> File "C:/Users/ddesantis/Dropbox/PythonScripts/CFD
> Models/ConversionModel.py", line 107, in
> res = Eq.sweep(dt=dt, solver=solver)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\term.py",
> line 236, in sweep
> solver = self._prepareLinearSystem(var=var, solver=solver,
> boundaryConditions=boundaryConditions, dt=dt)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\term.py",
> line 170, in _prepareLinearSystem
> buildExplicitIfOther=self._buildExplcitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\coupledBinaryTerm.py",
> line 122, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
> line 68, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
> line 68, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
> line 68, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\unaryTerm.py",
> line 99, in _buildAndAddMatrices
> diffusionGeomCoeff=diffusionGeomCoeff)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\abstractConvectionTerm.py",
> line 211, in _buildMatrix
> self.constraintB = -((1 - alpha) * var.arithmeticFaceValue *
> constraintMask * exteriorCoeff).divergence * mesh.cellVolumes
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
> line 1151, in __mul__
> return self._BinaryOperatorVariable(lambda a,b: a*b, other)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
> line 1116, in _BinaryOperatorVariable
> if not v.unit.isDimensionless() or len(v.shape) > 3:
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
> line 255, in _getUnit
> return self._extractUnit(self.value)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
> line 561, in _getValue
> value[..., mask] = numerix.array(constraint.value)[..., mask]
>
> IndexError: index 100 is out of bounds for axis 0 with size 100
>
> ___
> 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 ]
Re: casting implicit Boundary Conditions in FiPy
On Fri, Jul 1, 2016 at 2:44 PM, Gopalakrishnan, Krishnakumar <k.gopalakrishna...@imperial.ac.uk> wrote: > > Now, my question is, for the explicit term in the BC, can't I just define in > a regular Neumann BC style of definition ? , i.e. > phi.facegrad.constrain([alpha/ (1 - beta*dx/2)], mesh.facesRight]) . i.e. > will this ordinary Neumann correctly add up with the Implicit Boundary > condition term , implemented using the ImplicitSourceTerm method ? > > Or, should I treat this explicit BC term as part of the PDE itself, (like > the trick done for the implicit BC), ie. Do we need something along the lines > of fp.SourceTerm((mesh.faceNormals * explicitCoeff* > mesh.facesRight).divergence)) ? In either case, the diffusion > coefficient along the problematic boundary has been forced to be zero. I recommend trying both methods and implementing some 1D test cases to double check. I always do this no matter how much I trust the code. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
On Tue, Jun 28, 2016 at 12:53 PM, Abhilash Mathews <amath...@uwo.ca> wrote: > > eq1 = (TransientTerm(var=N1) == ImplicitSourceTerm(coeff=-2.*E1, var=P2) + > ImplicitSourceTerm(coeff= -2.*E2, var=P1) + ImplicitSourceTerm(coeff= > -1./(T1/TR), var=N1)) > > eq2 = (TransientTerm(var=P1) == ImplicitSourceTerm(coeff=2.*E2, var=N1) + > ImplicitSourceTerm(coeff= -1./(T2/TR), var=P1)) > > eq3 = (TransientTerm(var=P1) == ImplicitSourceTerm(coeff=2.*E1, var=N1) + > ImplicitSourceTerm(coeff= -1./(T2/TR), var=P2)) > > eq4 = (CentralDifferenceConvectionTerm(coeff = ones, var=E1) == > ImplicitSourceTerm(coeff=constant3, var=P2) + > ImplicitSourceTerm(coeff=-1./Ldiff, var=E1)) > > eq5 = (CentralDifferenceConvectionTerm(coeff = ones, var=E2) == > ImplicitSourceTerm(coeff=constant3, var=P1) + > ImplicitSourceTerm(coeff=-1./Ldiff, var=E2)) These are not great equations from FiPy's perspective as they are purely convective. However, if you want to persist with FiPy, you may want to use some sort of relaxation for the last two equations. One thing you could do is add in a TransientTerm and DiffusionTerm with small coefficients. Initially, just use coefficients of 1 to see if actually having those terms helps the numerical stability. If that helps, see if you can dial down those coefficients to get a more physical solution or use other techniques that maintain the correct physics. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
On Fri, Jun 17, 2016 at 5:36 PM, Gopalakrishnan, Krishnakumar <k.gopalakrishna...@imperial.ac.uk> wrote: > > My problem models the solid diffusion in a spherical particle. Matter > diffuses from the centre of the particle and reacts at the surface. This > is captured in a normalised 1-D domain with suitable equations and > co-ordinate scaling. Particle-centre is represented at x = 0 boundary, and > surface of the sphere represented by x=1 boundary. > > > > Now, my meshing algorithm is a special one. My inter-nodal distance (the > thickness of each sub-shell of the sphere) is such that, all of the internal > shells have equal volume. This is done so that mass is conserved within the > domain. This practise is stemming from my finite difference/finite element > colleagues who advocate this. Hi Krishna, I'm not sure what the thinking is of your colleagues, but the size of the elements has little or no impact on conservation. In finite volume, the equations are conservative to at least the precision of the linear solver (if not more) independently from the accuracy of the solution. Of course there could be an issue with source terms that I'm not seeing, but the diffusion and convection terms should be entirely conservative. > The problem is that, it is impossible to EXACTLY divide the shell into > integer number of iso-volume subshells. Thus, the scheme is chosen such > that, we get iso-volume shells for all inner sub-shells upto the last shell. > The last shell's thickness (dx) is obviously and purposefully made > ultra-small following the discussions and Dan's solution proposed in this > thread. Clearly, this last shell's volume differs from the rest of the inner > subshells. > > > > Am I violating any conservation laws here by doing it ? I know that finite > volume is a conservative method. But this question nevertheless nags me. I don't think so, what do the solutions show. Are the solutions conservative? > Is there any advantage to doing iso-volume subshells ( for the inner shells) > only to break this concept for the last shell ? I don't think so as long as the jump in the size of the cell volume is well within the precision of the solver. Three or four orders of magnitude should be okay. > Given that all the > solution dynamics happen at the surface (right boundary) , does a simple > geometric progression suffice with a very small dx at the right boundary > suffice ? Or is there any other optimal (structured) mesh generation > algorithm for choosing mesh sizes depending on problem-type, that I need to > refer to ? Don't know this one. Have you tried both methods and compared the solution? -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: bug in linearBicgstabSolver.py
Jamie, the new code is on the develop branch now https://github.com/usnistgov/fipy/blob/develop/fipy/solvers/scipy/linearBicgstabSolver.py It seemed to work for a small test case. Unfortunately, that solver isn't part of the test suite. On Wed, Jun 22, 2016 at 10:42 AM, James Pringle <jprin...@unh.edu> wrote: > thanks > Jamie -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: scipy's Delaunay output to fipy mesh2D object
James, this is awesome, thanks for posting. It would be a very good
idea to have a DelaunayMesh in FiPy that used Scipy for the
triangulation as it would give another way to make and test
unstructured meshes without the Gmsh dependency. Something like,
mesh = DelaunayMesh(points_in_2D)
In the long run, it would be good to change your code so that it
doesn'y loop in Python. I suspect that your code is quite inefficient
due to these loops. However, you have something that works, which is
great and it is only a one time cost at the beginning of the
simulation. Is efficiency an issue for you with this?
On Thu, Jun 16, 2016 at 4:29 PM, James Pringle <jprin...@unh.edu> wrote:
> For the use of others --
>
> Sometimes when dealing with arbitrary grids, it is useful to use the
> output of routines other than Gmsh to create the fipy mesh. In my case, I
> need to create a complex grid from ocean bathymetry, and I ended up using
> scipy.spatial.Delaunay to do so. With help from Daniel and Jonathan, I
> created a code to create a mesh2D object from the output of Delaunay. It is
> written to emphasize clarity over speed. I hope others find it useful.
>
> James Pringle
>
> #==
> #This code creates a fipy grid from the output of a Delaunay
> #triangulation. The code is written to prioratize clarity over speed.
>
> from pylab import *
> from numpy import *
> import fipy
> from scipy.spatial import Delaunay
>
>
> #make a simple example set of points to triangulate with delaunay
> points=array([[ 0., 0.],
>[ 0., 1.],
>[ 1., 0.],
>[ 1., 1.]])
> tri=Delaunay(points)
>
> #when this code is run as is, it should produce
> # faceVertexIDs=
> # [[3, 1, 0, 2, 0
> #[1, 0, 3, 3, 2]]
> #and
> # cellFaceIds=
> #[[0, 3],
> # [1, 2],
> # [2, 4]]
>
>
> #now create the arrays that fipy.mesh.mesh2D needs. vertexCoords is
> #of shape (2,N), where N is the number of points, and contains the
> #coordinates of the vertices. It is equivalent to tri.points in
> #content, but reshaped.
> print 'Making vertexCoords'
> vertexCoords=tri.points.T
>
> #faceVertexID is of shape (2,M) where M is the number of faces/edges
> #(faces is the fipy terminology, edges is the scipy.spatial.Delaunay
> #terminology). faceVertexID contains the points (as indices into
> #vertexCoords) which make up each face. This data can be extracted
> #from simplices from Delaunay, which contains the faces of each
> #triangle. HOWEVER, simplices contains many repeated faces, since each
> #triangle will border on others. Only one copy of each face should be
> #inserted into faceVertexID.
>
> print 'Making faceVertexIDs'
> faceIn={} #this is a dictionary of all faces that have already been
> #inserted into the faceVertexID, with the key a tuple of
> #indices, sorted
> faceVertexIDs_list=[] #structure to build
>
> for ntri in xrange(tri.simplices.shape[0]):
> for nface in xrange(3):
> if nface==0:
> face=tri.simplices[ntri,[0,1]]
> elif nface==1:
> face=tri.simplices[ntri,[1,2]]
> else:
> face=tri.simplices[ntri,[2,0]]
> faceKey=tuple(sort(face))
> if not (faceKey in faceIn):
> faceIn[faceKey]=True
> faceVertexIDs_list.append(face)
>
> faceVertexIDs=array(faceVertexIDs_list).T
>
> #now create cellFaceIDs of shape (3,P) where P is the number of cells
> #in the domain. It contains the faces (as indices into faceVertexIDs)
> #that make up each cell.
>
> #first create dictionary with a key made up of a sorted tuple of vertex
> #indices which maps from a face to its location in faceVertexIDs
> print 'Making cellFaceIDs'
> faceMap={}
> for nface in xrange(faceVertexIDs.shape[1]):
> faceKey=tuple(sort(faceVertexIDs[:,nface]))
> faceMap[faceKey]=nface
>
> #now loop over simplices, find each edge/face, and map from points to
> #faces
> cellFaceIDs=tri.simplices.T*0
> for ntri in xrange(tri.simplices.shape[0]):
> for nface in xrange(3):
> if nface==0:
> face=tri.simplices[ntri,[0,1]]
> elif nface==1:
> face=tri.simplices[ntri,[1,2]]
> else:
> face=tri.simplices[ntri,[2,0]]
> faceKey=tuple(sort(face))
> cellFaceIDs[nface,ntri]=faceMap[faceKey]
>
> #ok, now instantiate a mesh2D object with this data
> print 'Making mesh'
> mesh=fipy.meshes.mesh2D.Mesh2D(vertexCoords,faceVertexIDs,cellFaceIDs)
>
>
>
> ___
> 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 ]
Re: how to create large grid with holes and odd geometry
I think that you need to instatiate one of these https://github.com/usnistgov/fipy/blob/develop/fipy/meshes/mesh2D.py#L70 with the correct arrays. For examples, >>> points = [[0,0], [0, 1], [1, 0], [1, 1]] >>> from scipy.spatial import Delaunay >>> tri = Delaunay(points) >>> tri.points array([[ 0., 0.], [ 0., 1.], [ 1., 0.], [ 1., 1.]]) >>> tri.simplices array([[3, 1, 0], [2, 3, 0]], dtype=int32) So the vertices are just the points but reshaped to (2, N) from (N, 2). The face to vertex array would be [[3, 1, 0, 2, 0 [1, 0, 3, 3, 2]] Its all the faces or edges in 2D, and the cell to face ID would be [[0, 3], [1, 2], [2, 4]] Then you can instatiate the Mesh2D object with those three arrrays. Look at https://github.com/usnistgov/fipy/blob/develop/fipy/meshes/mesh2D.py#L269 to see how to do it as well. Notice that the test example has both quads and triangles (the -1 values are masked). On Wed, Jun 15, 2016 at 1:29 PM, James Pringle <jprin...@unh.edu> wrote: > Well, I am motivated to give it a go, since I only have the summer to make > progress on project and it is blocking my research progress right now. Can > you give me a pointer to where the appropriate quantities are defined? I can > certainly write code to make the transformations, but it is a bit hard > without understanding precisely what is defined in the Mesh2D object. I have > made a simple Mesh2D object, but I am not sure which of the attributes, etc, > are absolutely required, or how they are precisely defined. > > Perhaps most useful for me would be > > a definition of which parts of the Mesh2D object must exist > and the format of those parts, in particular the a face to vertex array and > a > cell to face array > > Or you could point me to the appropriate part of the Gmsh code so I can go > from there. I presume poking around in fipy.Gmsh2D would be a good place to > start? Is there a better place to start? > > I would love any documentation on the details of the Mesh2D object. > > Thanks, > Jamie Pringle > University of New Hampshire > > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler <daniel.wheel...@gmail.com> > wrote: >> >> Hi Jamie, >> >> There is no automated way to make a FiPy mesh from Scipy's Delaunay >> object. The problem is that FiPy needs a face to vertex array and a >> cell to face array while the Delaunay object just has the cell to >> vertex array. The functionality to extract the face to vertex array >> and cell to face array is in FiPy because it must be done for Gmsh, >> however, it is not abstracted in so that it can be reused. >> >> It is certainly possible to make the correct arrays with Numpy and >> pass them to the Mesh2D object, but it's a bit of work to write the >> code. If I find some time I might give it a go, but I don't know when >> I will get to it. >> >> Cheers, >> >> Daniel >> >> On Tue, Jun 14, 2016 at 4:15 PM, James Pringle <jprin...@unh.edu> wrote: >> > Daniel et al. -- >> > >> > As referenced earlier in this thread, I have a complex domain with >> > holes >> > in it; I now have it broken up into triangles, based on the Delaunay >> > package >> > in SciPy. I have >> > >> > Locations of all vertex coordinates, >> > list of vertices that make up faces >> > list of faces that make cells >> > list of faces that make up all the internal and external boundaries. >> > >> > Can you point me to any code or documentation that would help me >> > understand >> > how to make this into a Mesh2D object? I am having a devil of a time >> > figuring out from the manual online. The best would be something that >> > used >> > the output of either the triangles or scipy.spatial.Delaunay() packaged. >> > My >> > equation is of the form >> > >> > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) >> > >> > >> > and I can get the coefficients A(x,y) and B(x,y) on either the faces or >> > in >> > the cell centers are needed. >> > >> > Thank you. >> > Jamie Pringle >> > University of New Hampshire >> >> -- >> Daniel Wheeler >> ___ >> fipy mailing list >> fipy@nist.gov >> >> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy=CwICAg=c6MrceVCY5m5A_KAUkrdoA=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclilpum_QCUOE= >> [ NIST internal ONL
Re: casting implicit Boundary Conditions in FiPy
On Wed, Jun 15, 2016 at 7:27 AM, Gopalakrishnan, Krishnakumar <k.gopalakrishna...@imperial.ac.uk> wrote: > > Dan. I was able validate that your code correctly implements the Implicit > Neumann Boundary Condition (in 1-D). > > Here is a link to the plot of the solution after a transient simulation for > 0.2 seconds. The plot on the left is generated in Mathematica, and that to > the right is that generated by the FiPy’s matplotlib viewer class. The two > are identical. > > https://imperialcollegelondon.box.com/s/lb8v1iqxb3rqhxsphn9cmhwcmh5jzpiz Awesome that it worked. I think that you're right though regarding first versus second order boundary condition. I think it is first order in space. -- Daniel Wheeler ___ 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 create large grid with holes and odd geometry
Hi Jamie, There is no automated way to make a FiPy mesh from Scipy's Delaunay object. The problem is that FiPy needs a face to vertex array and a cell to face array while the Delaunay object just has the cell to vertex array. The functionality to extract the face to vertex array and cell to face array is in FiPy because it must be done for Gmsh, however, it is not abstracted in so that it can be reused. It is certainly possible to make the correct arrays with Numpy and pass them to the Mesh2D object, but it's a bit of work to write the code. If I find some time I might give it a go, but I don't know when I will get to it. Cheers, Daniel On Tue, Jun 14, 2016 at 4:15 PM, James Pringle <jprin...@unh.edu> wrote: > Daniel et al. -- > > As referenced earlier in this thread, I have a complex domain with holes > in it; I now have it broken up into triangles, based on the Delaunay package > in SciPy. I have > > Locations of all vertex coordinates, > list of vertices that make up faces > list of faces that make cells > list of faces that make up all the internal and external boundaries. > > Can you point me to any code or documentation that would help me understand > how to make this into a Mesh2D object? I am having a devil of a time > figuring out from the manual online. The best would be something that used > the output of either the triangles or scipy.spatial.Delaunay() packaged. My > equation is of the form > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > > and I can get the coefficients A(x,y) and B(x,y) on either the faces or in > the cell centers are needed. > > Thank you. > Jamie Pringle > University of New Hampshire -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
On Sun, Jun 12, 2016 at 8:14 AM, Gopalakrishnan, Krishnakumar <k.gopalakrishna...@imperial.ac.uk> wrote: > > 1. The backward Euler method is used in formulating the equation, > "n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) " . But Backward Euler > is a first order method, isn't it ? I am a bit confused about the second > order accuracy statement. I'm confused by your question. Backward Euler refers to the time stepping scheme, which is indeed first order. The spatial discretization scheme is second order and using the boundary condition trick that I demonstrated is second order and also implicit. > 2. In, the equation, "eq = (fp.TransientTerm() == > fp.DiffusionTerm(diffCoeff) - \ > > fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * > mesh.facesRight).divergence)) " , > > from what I understand the 'ImplicitCoeff' is zeroed out in all other > locations other than the right boundary by the boolean method > 'mesh.facesRight' , right ? Exactly, mesh.facesRight is just evaluates to a boolean array. There is also mesh.exteriorFaces which gets you all exterior faces and then you can also do logical operations based on physical location to get more specific masks to capture boundaries or internal regions. > If this is being implemented, do we have the restriction to use fixed dx ? I guess not if you're careful, but it could be a bit fiddly getting the correct distance. > Although dx appears in the Implicit term, it can just be set to the value > of 'dx' of the last node (i.e. length of the last segment in a > variable-mesh sized 1D geometry file).Am I missing something here ? No, but you need to know how to get the correct distance. There are arrays of cell to face distances, but I'd have to think about how to get that working right. > 3. Also, I am a bit confused about the negative signs used in the equation > above.The implicitCoeff is given a negative sign, and the equation's > ImplicitSourceTerm is also given a negative sign. From what I undertand, > this is what we intend to do ? I was just flipping signs to make it work so I'm not sure. Maybe it would be better to have them both positive and it would make more sense that way. > Assign the diffusion coefficient D throughout the domain, using the > faceVariable method. > Manually set it to zero at the right boundary > Then **add** the Diffusion achieved by the div.(D grad( phi)) by using the > ImplicitSourceTerm method, and ensure that this effect gets added only at > the right boundary. You got it. > If this is the workflow, then perhaps we just have to declare the > ImplicitSourceTerm to be positive, and use the addition sign instead of > subtraction sign for the implicitSourceTerm in the equation ? Yup, I think so. > 5. Taking one massive step to get to the answer. Here is a big question > (due to my poor knowledge in the subject) > > I understand that the implicitness enables us to take these large time-steps > due to the affine and A-stable properties of Implicit Methods.But are we > required to do this ? Oh no. You can still take as small a time step as necessary to capture the dynamics at the time scale of interest. It's just that you're not limited by the horrible dx**2 / D time step limit associated with explicit diffusion. Typically it's best to resolve to some dx based time scale (rather than dx**2) that captures an interface being tracked or some feature of interest moving around. >My problem is that, I am solving a coupled system > of PDEs and my other equations, require me to step very very slowly, due to > the inherent stiffness of my system. The PDEs describe processes of varying > time-constants. So, does taking the timestep to be small inhibit us > in any way ? The only down side to small time steps is that their small :-) > Once again, thanks a lot for your solution. Based on the discusssions with > Ray and from this particular method from your posting to the list, I > think I shall be able to solve the problem in 1D. The issues you pointed > out about using this approach in 2D went right over my head. I have to > think about this, when I implement them perhaps in a week or so. Good luck with 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 ]
Re: casting implicit Boundary Conditions in FiPy
>> setting phi.faceGrad.constrain([k*phi], mesh.facesRight), and see if it
>> will work.
>>
>>
>>
>> Thanks for pointing this out.
>>
>>
>>
>>
>>
>> Krishna
>>
>>
>>
>> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov
>> <fipy-boun...@nist.gov>] *On Behalf Of *Raymond Smith
>> *Sent:* 08 June 2016 23:36
>> *To:* fipy@nist.gov
>> *Subject:* Re: casting implicit Boundary Conditions in FiPy
>>
>>
>>
>> Hi, Krishna.
>>
>> Just to make sure, do you mean that the boundary condition is a
>> derivative with respect to the spatial variable or with respect to time
>> as-written? If you mean spatial, such that d\phi/dx = k*phi, have you tried
>>
>> phi.faceGrad.constrain(k*phi) and that didn't work?
>>
>> If you mean that its value is prescribed by its rate of change, then I'm
>> not sure the best way to do it. Could you maybe do it explicitly?
>> - Store the values from the last time step with hasOld set to True in
>> the creation of the cell variable
>> - In each time step, calculate the backward-Euler time derivative
>> manually and then set the value of phi with the phi.constrain method.
>>
>>
>>
>> Ray
>>
>>
>>
>> On Wed, Jun 8, 2016 at 6:26 PM, Gopalakrishnan, Krishnakumar <
>> k.gopalakrishna...@imperial.ac.uk> wrote:
>>
>> I am trying to solve the standard fickean diffusion equation on a 1D
>> uniform mesh in (0,1)
>>
>>
>>
>> $$\frac{\partial \phi}{\partial t} = \nabla.(D \nabla\phi)$$
>>
>>
>>
>> with a suitable initial value for $\phi(x,t)$.
>>
>>
>>
>> The problem is that, one of my boundary conditions is *implicit*, i.e.
>> is a function of the field variable being solved for.
>>
>>
>>
>> $ \frac{\partial\phi}{\partial t} = k \phi $ , at the right boundary
>> edge, k = constant
>>
>>
>>
>> The left BC is not a problem, it is just a standard no-flux BC.
>>
>>
>>
>> How do I cast this *implicit BC* in FiPy ? Any help/pointers will be
>> much appreciated.
>>
>>
>>
>>
>>
>> Best regards
>>
>>
>>
>> Krishna
>>
>> Imperial College London
>>
>>
>> ___
>> 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 ]
>>
>>
>>
>> ___
>> 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 ]
>
>
--
Daniel Wheeler
___
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 create large grid with holes and odd geometry
Hi Jamie, The simplest way might be just to take a basic grid and just "switch off" the simulation outside of the black zone. This can be achieved by including a transient term in the equation and setting the coefficient to be large in the white zone and zero or something small in the black zone. That would be an easy first attempt to solve the problem. Of course this approach has a lot of redundancy, but would be a good first step. You may have to interpolate the A and B values to a coarser grid for this approach to be feasible. It is possible in FiPy to combine grid's into odd shaped non-rectangular grids. You could do this with say 20 or 30 smaller grids to better isolate the black region. I still think that you want to interpolate the values of A and B so you have control over the grid density. You don't necessarily need a grid density at the sample density. Also, you will still need the trick with the Transient Term. Furthermore, the grids must align and not overlap. It will take quite a bit of programming to set this up so that you have control over the grid density and the number of subgrids to isolate the black region. Gmsh will probably be happy meshing the black region, but will probably use triangles rather than a grid. To summarize use the first approach to make things work right and then refine with the second approach. Hope it helps. On Tue, Jun 7, 2016 at 1:46 PM, James Pringle <jprin...@unh.edu> wrote: > Dear mailing list & developers -- > > I am looking for hints on the best way to proceed in creating a > grid/mesh for a rather complex geometry. I am just looking for which method > (Gmsh or something else?) to start with, so I can most efficiently start > coding without exploring blind alleys. > > I am solving an elliptic/advective problem of the form > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > where Psi is the variable to solve for, and A(x,y) and B(x,y) are > coefficients known on a set of discrete points shown as black in > https://dl.dropboxusercontent.com/u/382250/Grid01.png . The black appears > solid because the grid is dense. > > The locations of the points where the coefficients are known define the > grid. The number of points is large (911130 points) and they are evenly > spaced where they exist. Note that there are holes in the domain that > represent actual islands in the ocean. > > I am happy to keep the resolution of the grid/mesh equal to the spacing of > the points where the coefficients are known. > > What is the best way to approach creating a grid for this problem? I would > love code, of course, but would be very happy with suggestions of the best > way to start. > > Thanks > Jamie Pringle > University of New Hampshire > > ___ > 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 ]
Re: transient term with coefficient
Hi Francisco, In FiPy, it is necessary to add both a TransientTerm and a source term for terms such as, density (times) (\partial f[x]/\partial t) So, the code will look something like this, TransientTerm(var=var, coeff=density) - field * (density - density.old) / dt There may be ways to improve the convergence with implicit sources, but I haven't investigated that. For example, the "field * density.old / dt" part could become implicit in "field" and the the "field * density / dt" could be implicit for "density" and coupled with an equation solving for density. Cheers, Daniel On Mon, May 16, 2016 at 1:30 PM, Francisco Vega Reyes <fv...@unex.es> wrote: > Hello, > > a very common term in the Navier-Stokes equations for a gas is of the form: > > density (times) (\partial f[x]/\partial t), > > where f[x] is a hydrodynamic field (for instance, temperature). however, it > appears in the Fipy manual that the transient terms that express time > derivatives have no multiplying coefficients. I’d like to solve the transient > Navier Stokes equations for a gas without having to divide the whole balance > equations by density, which is of course also an unknown. > > How can I express then a term of the form above in Fipy? > > Thank you very much, > > > Dr. Francisco Vega Reyes > > Departamento de Física, > Universidad de Extremadura, > 06071 Badajoz, Spain > > fv...@unex.es > > > > > > ___ > 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 ]
Re: IndexError: index is out of bounds for axis 1 with size
One possible way to handle this would be to have Pu and Pd defined
everywhere, but only have meaning in the edge cells. The coefficients C3
and C4 would only have a non-zero value in the end cells. Move the end
cells so the cell centers are on the boundaries. The boundary conditions
for the Pf and Pk equations can then be source terms in the edge cells that
constrain the edge values. The equations can then be tightly coupled
implicitly. I think this will all work quite well. The important thing to
know is that a single cell value can be constrained with the application of
a large implicit and explicit source, e.g.,
x, y = mesh.cellCenters
big_coeff = CellVariable(mesh=mesh, value=0.)
big_coeff[x < dx] = 1e+20
big_coeff[x > L - dx] = 1e+20
TransientTerm() == big_coeff * value - ImplicitSourceTerm(big_coeff)
if `L` is the length of the mesh and `dx` is the grid spacing. The edge
cells will be constrained to `value`. `value` can be one of the
CellVariables that are in the other equations.
I hope that helps you get started.
On Mon, May 9, 2016 at 11:57 AM, Mohammad Kazemi <kazem...@gmail.com> wrote:
> I forgot to mention that the [image: P_f] and [image: P_k] values in eqpd
> and eqpu are actually as follow,
>
> eqpd:
>
> [image: -\frac{\partial P_d}{\partial t} = C_3 (P_f (0,t) - P_d(t)) + C_4
> (P_k(0,t)-P_d(t))]
>
> and eqpu:
>
> [image: \frac{\partial P_u}{\partial t} = C_3 (P_f (L,t) - P_u(t)) + C_4
> (P_k(L,t)-P_u(t))]
>
>
> On Mon, May 9, 2016 at 11:19 AM, Mohammad Kazemi <kazem...@gmail.com>
> wrote:
>
>> Thanks for your response. Here are my equations:
>>
>> eqpk:
>>
>> [image: \frac{\partial P_k}{\partial t} = D \frac{\partial^2
>> P_k}{\partial x^2}]
>> with IC and BCs:
>> [image: P_k (x,0) = Pi]
>> [image: P_k(0,t) = P_d(t)]
>> [image: P_k(L,t) = P_u(t)]
>>
>> eqpf:
>>
>> ∂Pf/∂t = C1 ∂2 Pf/∂x2 - C2 (Pf-Pk)
>>
>> with IC and BCs:
>>
>> [image: P_f (x,0) = Pi]
>> [image: P_f(0,t) = P_d(t)]
>> [image: P_f(L,t) = P_u(t)]
>>
>> eqpd:
>>
>> [image: - \frac{\partial P_d}{\partial t} = C_3 (P_f - P_k)+C_4 (P_k-P_d)]
>>
>> with IC:
>>
>> [image: P_d (0) = Pi]
>>
>> eqpu:
>>
>> [image: \frac{\partial P_u}{\partial t} = C_3 (P_f - P_u)+C_4 (P_k-P_u)]
>>
>> with IC:
>>
>> [image: P_u (0) = Pi]
>>
>> So [image: P_k] and [image: P_f] are functions of space and time while
>> [image:
>> P_d] and [image: P_u] are only function of time. How should i define a
>> variable which is only a function of time? If I want to have schematic of
>> problem, it is as below.
>>
>>
>>
>>
>>
>> The pressure in tanks [image: P_u] and [image: P_d] are uniform and only
>> changes with time.
>>
>> I appreciate it,
>>
>>
>>
>>
>> On Mon, May 9, 2016 at 9:43 AM, Guyer, Jonathan E. Dr. (Fed) <
>> jonathan.gu...@nist.gov> wrote:
>>
>>> By declaring Pd and Pu as CellVariables, you are declaring them to be
>>> functions of space. That's what a CellVariable is.
>>>
>>> Furthermore, eqpd and eqpu depend on Pf and Pk, so that necessitates
>>> that Pd and Pu are functions of space.
>>>
>>> Can you show us the math for the equations you are trying to solve?
>>>
>>> > On May 7, 2016, at 1:12 PM, Mohammad Kazemi <kazem...@gmail.com>
>>> wrote:
>>> >
>>> > Thank you. It is running now, but I have a question. Pf and Pk are
>>> functions of location (x) and time (t). However, Pd and Pu are only
>>> functions of time (they represent pressure at downstream and upstream
>>> tanks). When I solve this problem, I will get a vector for Pd and Pu (at
>>> last timestep) which varies with location. For example, for Pu,
>>> >
>>> > print Pu
>>> >
>>> > [ 7171802.05029459 7171719.44929532 7171632.68237858
>>> 7171538.95796536
>>> > 7171435.46595385 7171319.37283423 7171187.81897787
>>> 7171037.91833804
>>> > 7170866.76056175 7170671.41527712]
>>> >
>>> > Which shows the Pu values at 10 different mesh blocks. Is it related
>>> to how I define my cell variables?
>>> >
>>> > I appreciate your helps.
>>> >
>>> > Bests,
>>> >
>>> > On Fri, May 6, 2016 at 5:21 PM, Daniel Wheeler <
>>> daniel.wheel...@gmail.com> wrote:
>>> > The following runs for me with a few changes.
>>> >
>>> > On Fri, May 6, 2016 at 1:25 PM, Moham
Re: IndexError: index is out of bounds for axis 1 with size
The following runs for me with a few changes. On Fri, May 6, 2016 at 1:25 PM, Mohammad Kazemi <kazem...@gmail.com> wrote: > > Pk.constrain([Pd], mesh.facesLeft) > Pk.constrain([Pu], mesh.facesRight) Pk.constrain(Pd.faceValue, mesh.facesLeft) Pk.constrain(Pu.faceValue, mesh.facesRight) > Pf.constrain([Pd], mesh.facesLeft) > Pf.constrain([Pu], mesh.facesRight) Pf.constrain(Pd.faceValue, mesh.facesLeft) Pf.constrain(Pu.faceValue, mesh.facesRight) This makes sense since the Pf and Pk are scalars and [Pd] and [Pu] are effectively vectors. Also the face value is required since the constraint is on the faces. I hope this helps. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: IndexError: index is out of bounds for axis 1 with size
On Fri, May 6, 2016 at 12:09 PM, Daniel Wheeler <daniel.wheel...@gmail.com> wrote: > Hi Mo, Apologies for that. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: IndexError: index is out of bounds for axis 1 with size
Hi Mo, I tried running your code but I got "name 'Pf' is not defined" -- a different error. On Fri, May 6, 2016 at 12:41 AM, Mohammad Kazemi <kazem...@gmail.com> wrote: > > CODE: > > from __future__ import division > from numpy import * > from pylab import * > from matplotlib import rc > from fipy import * > > INITIALLIZATION > ### > L = 0.0508 #+0.0021+0.0014 > nx = 10.0 > dx = L / nx > Pdi=6894757.29 > Pui=7170547.58 > D=1e-6 > tmax = 1000 > sigma = 400 > k0 = 10e-16 > T = 300 > M = 0.0399 > R = 8.314 > muf=1.96e-5 > hf=30e-6 > phif=0.25 > A1=1.14 > A2=0.28 > wkf=4 > wkd=4 > wku=4 > Kn=sqrt((3.14*R*T/(2*M))*muf/(Pf*hf)) At this point "Pf" is not defined. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: understanding convection terms
On Tue, Apr 26, 2016 at 10:57 AM, Kris Kuhlman <kristopher.kuhl...@gmail.com> wrote: > Daniel, > > Thank you. I am a bit surprised that the CentralDifference basically matches > the hybrid method, and is more accurate than upwind. Remember that central difference is second order accurate. The other schemes are first order. Of course it's unstable when the Peclet number is above 2 that's why we use the other schemes. Basically, the schemes (hybrid, exponential etc) are ways to swap between central difference and upwind based on the local Peclet number. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: understanding convection terms
Hi Kris,
Good to hear from you again and a very nicely coded example.
I think it's a solver issue. I tried the LU solver and it gives
perfect results. See
https://gist.github.com/wd15/affe4d82cc2a189d894a7d774e4bc00b.
This might suggest that we need to change the default solver when
using the central difference scheme.
Cheers,
Daniel
On Mon, Apr 25, 2016 at 11:54 AM, Kris Kuhlman
<kristopher.kuhl...@gmail.com> wrote:
> I am trying to understand the convection terms available in fipy through a
> simple steady-state problem. I am surprised at the divergence of the
> CentralDifferenceConvectionTerm, is this to be expected? As the
> discretization in the mesh is made finer, the solution gets worse!?
>
> The problem is: \frac{\partial^2 u}{\partial x^2} - v \frac{\partial
> u}{\partial x} = 0
>
> The script at the link below compares the solution using different
> ConvectionTerms, and plots the figures attached for different values of nx.
> For larger values of v, the solution diverges even more and becomes
> oscillatory.
>
> https://gist.github.com/klkuhlm/07f9eaf52b24e103f60ae213c0944c21
>
> Is this expected behavior?
>
> Kris
>
>
> ___
> 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 ]
Re: Thermal Diffusion with Source
On Tue, Apr 19, 2016 at 11:10 PM, John Leeman <kd5...@gmail.com> wrote: > > I was just worried I didn’t grasp the way I needed to interact with FiPy. An > example along this vein could be a useful addition in the docs possibly? John, an example would be very welcome. If you would like to just make it an IPython notebook (rather than doctest/restructured text script) then please do that and we'll try and link to it from the documentation. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: squared partial derivatives
On Tue, Apr 19, 2016 at 11:22 AM, Francisco Vega Reyes <fv...@unex.es> wrote: > Hello, > > I am trying to solve a compressible flow problem (Navier Stokes > equations for a gas). > > Can a term like: > > (\partial(f(x,y))/\partial y) * (\partial(f(x,y))/\partial y) > > be represented in FiPy ? Hi Francisco, You can represent any term explicitly in FiPy, including (f_y)^2. However, that sort of term might arise from a conservation law of the form (f f_y)_y. If that is the case, it is good to try and implement terms in their conservative form. In the case of (f f_y)_y, "f" will be a coefficient in a diffusion term. Cheers, Daniel -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: fipy diffusion setup
On Tue, Apr 19, 2016 at 4:10 PM, Phil Battle <bat...@advr-inc.com> wrote: > Hi Daniel > Thank you for following up- > I did try your suggestion, ( C is my concentration dependent Cell variable) > > Dx=Dpe_x*(a_x+((1.-a_x)/(b_x*C+g_x))) > Dz=Dpe_z*(a_z+((1.-a_z)/(b_z*C+g_z))) > diff = Dx.dot([[1,0],[0,0]])+Dz.dot([[0,0],[0,1]]) > > and program runs and "looks correct" so good to go!! Great news. > But before your suggestion I was trying to create the suggested rank 2 > diffusion coefficient a different way and I did not get far before hitting a > road block- If you have time to answer this specific question it may help > with understanding the structure of variables in this program > > I created a 2x2 Cell Variable on a 2x2 2Dim mesh. I can show that code if > necessary, but the result for the initialized CellVariable is shown below > > print(D.value) > > [[[ 0. 2. 4. 6.] > [ 1. 1. 1. 1.]] > > [[ 1. 1. 1. 1.] > [ 1. 1. 1. 1.]]] > > > Now the part that I am having trouble with. How can I change the 6 to a 0.? I think you just need to do D[0, 0, 3] = 6. However, it isn't always a good practice to access the spatial index since if the grid changes then the index has a different spatial location > The closet I can come is changing all the values at that mesh position to 0. > that is using > > D.setValue(0,where=[[0,0,0,1]]) I'm not sure why the [[0, 0, 0, 1]] does what it does. Numpy does the following: ~~~ In [9]: a = np.arange(16).reshape([2, 2, 4]) In [10]: a Out[10]: array([[[ 0, 1, 2, 3], [ 4, 5, 6, 7]], [[ 8, 9, 10, 11], [12, 13, 14, 15]]]) In [12]: a[[[0, 0, 0, 1]]] = -1 In [13]: print a [[[-1 -1 -1 -1] [-1 -1 -1 -1]] [[-1 -1 -1 -1] [-1 -1 -1 -1]]] ~~~ If you're trying to do fancy indexing with FiPy, you might want check with Numpy first so that you can trust the FiPy result. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: fipy diffusion setup
Phil, just to follow up and complete the answer. I neglected to explain how to create a rank 2 CellVariable from two rank 0 CellVariables to preserve the dependency. The dot operator helps with that. Given two rank 0 CellVariables, v00 and v11, then to get the rank 2 CellVariable use, vv = v00.dot([[1, 0], [0, 0]]) + v11.dot([[0, 0], [0, 1]]) On Mon, Apr 18, 2016 at 3:33 PM, Daniel Wheeler <daniel.wheel...@gmail.com> wrote: > On Mon, Apr 18, 2016 at 12:33 PM, Phil Battle <bat...@advr-inc.com> wrote: >> Ok, below works as I would expect ( note Dx and Dz are constants) > > Thanks for that. I can confirm that it's broken in Python 2.7 as well > when the diffusion coefficient is [[Dx, Dy]]. > > > -- > Daniel Wheeler -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
