Re: Solving advection, segregation and diffusion equation
> On Sep 14, 2016, at 12:49 PM, Zhekai Deng > wrote: > > So, in my previous one, I only have constraints on left (fixed 0.5 inlet > concentration), top and bottom ( the boundary conditions that I applied), but > not boundary condition on the right. I have implemented an additional one on > the right to account for the outlet: > > phi.faceGrad.constrain(velocityVector*phi.faceValue,where = mesh.facesRight) > > but it seems does not help very much. I was surprised to see that several "corrections" I applied to your script didn't seem to change things very much. This problem is strangely insensitive to things that I think should matter. In contrast, moving the definition of eq out of the sweep loop made a huge difference and I don't understand why. > I am not sure what does "velocity is backwards someplace" mean. From the way > I defined it, the $u \hat{i} = \tilde{z}$ and $u > \hat{j} = 0$. Both are independent of solution variable. Would it be possible > to clarify what does "velocity is backwards" mean? This was just idle speculation that maybe a velocity or shear points left when it should point right (up vs. down), etc. and that maybe that's what causes the rippling. > After the discussion with my colleague, we think it is better to understand > this as a vector term that only act on the normal direction ($\hat{j}$) of > the flow domain. The reason for this is that the shear rate term is gradient > of velocity field which is a tensor. However, we are not entire sure the > functional dependence of segregation velocity on shear rate tensor. That's > why we usually call it as local shear rate (du/dz \hat{j} ), and I should > have been clear on this issue in the beginning. Thanks for clarifying. > In summary, it does seems that the initial sweeps give the right "trend" > toward the solution. However, I am still not sure about why it is not > convergent. Any suggestion on how to proceed to the next step ? I have > attached the newest version of the code. 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. - Jon ___ 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
Thank you for the reply. I try to answer any confusion in below. 1. I think what you have is fine. Explicitly, I'd write: phi.faceGrad.constrain(..., where=mesh.facesTop | mesh.facesBottom) See http://www.ctcms.nist.gov/fipy/documentation/USAGE.html# applying-spatially-varying-boundary-conditions Thanks. It has been implemented in this way now. 2. For inlet/outlet conditions, please see the comment at http://www.ctcms.nist.gov/fipy/documentation/USAGE.html# applying-outlet-or-inlet-boundary-conditions If you impose constraints, then FiPy is inlet/outlet, not no-flux. So, in my previous one, I only have constraints on left (fixed 0.5 inlet concentration), top and bottom ( the boundary conditions that I applied), but not boundary condition on the right. I have implemented an additional one on the right to account for the outlet: phi.faceGrad.constrain(velocityVector*phi.faceValue,where = mesh.facesRight) but it seems does not help very much. I've made some changes to your script that result in the initial sweeps looking much more like your matlab solution, but it is still not convergent. I wonder at this point if a velocity is backwards someplace. Thank you for the change. Well, even after I added the right side outlet boundary condition, the sweep solution still looks similar to what you have. In the first initial sweeps, it looks like the matlab solution, but it is not convergent. However, the residual seems does not explore as quickly as before by adding the outlet condition. I am not sure what does "velocity is backwards someplace" mean. From the way I defined it, the $u \hat{i} = \tilde{z}$ and $u \hat{j} = 0$. Both are independent of solution variable. Would it be possible to clarify what does "velocity is backwards" mean? I don't understand why $\tilde{\dot{\gamma}} = 1 \hat{\j}$. Shouldn't it be $\tilde{\dot{\gamma}} = 1 \hat{\i}$, given that $\tilde{u} = \tilde{z} \hat{\i}$? Thanks for bring it up. After the discussion with my colleague, we think it is better to understand this as a vector term that only act on the normal direction ($\hat{j}$) of the flow domain. The reason for this is that the shear rate term is gradient of velocity field which is a tensor. However, we are not entire sure the functional dependence of segregation velocity on shear rate tensor. That's why we usually call it as local shear rate (du/dz \hat{j} ), and I should have been clear on this issue in the beginning. Is, perhaps, the definition $\tilde{\dot{\gamma}} \equiv \frac{\partial \tilde{u}_x}{\partial \tilde{z}} \hat{\j}$ where $\tilde{u} \equiv \tilde{u}_x \hat{\i} + \tilde{u}_y \hat{\j}$ Yes, this is the definition we are currently following. In summary, it does seems that the initial sweeps give the right "trend" toward the solution. However, I am still not sure about why it is not convergent. Any suggestion on how to proceed to the next step ? I have attached the newest version of the code. Best, Zhekai On Tue, Sep 13, 2016 at 6:35 PM, Guyer, Jonathan E. Dr. (Fed) < jonathan.gu...@nist.gov> wrote: > Zhekai - > > 1. I think what you have is fine. Explicitly, I'd write: > > phi.faceGrad.constrain(..., where=mesh.facesTop | mesh.facesBottom) > > See http://www.ctcms.nist.gov/fipy/documentation/USAGE.html# > applying-spatially-varying-boundary-conditions > > > 2. For inlet/outlet conditions, please see the comment at > > http://www.ctcms.nist.gov/fipy/documentation/USAGE.html# > applying-outlet-or-inlet-boundary-conditions > > If you impose constraints, then FiPy is inlet/outlet, not no-flux. > > > I've made some changes to your script that result in the initial sweeps > looking much more like your matlab solution, but it is still not > convergent. I wonder at this point if a velocity is backwards someplace. > > I don't understand why $\tilde{\dot{\gamma}} = 1 \hat{\j}$. > Shouldn't it be $\tilde{\dot{\gamma}} = 1 \hat{\i}$, given that $\tilde{u} > = \tilde{z} \hat{\i}$? > > Is, perhaps, the definition > $\tilde{\dot{\gamma}} \equiv \frac{\partial \tilde{u}_x}{\partial > \tilde{z}} \hat{\j}$ > where > $\tilde{u} \equiv \tilde{u}_x \hat{\i} + \tilde{u}_y \hat{\j}$ > > (I *really* strongly recommend not writing vectors like scalars) > > - Jon > > > On Sep 13, 2016, at 5:45 PM, Zhekai Deng northwestern.edu> wrote: > > > > Hello All, > > > > I am trying to use Fipy to solve advection, segregation and diffusion > equation in 2D. I have attached the pdf to explain in details about my > equation, and corresponding fipy code that I have problem with it. > > > > The solution I got has "kind of similar shape" to the same equation that > I solved using MATLAB(the code is also attached, I have verified that > MATLAB is solving correctly). However, the code I have in Fipy does not > converge, and the residual jump between large and small number. > > > > Also,I have a few things that I am not entire sure. > > > > 1. In 2D, how could we specify the which direction of the > "faceGrad.constrain" t
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) 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 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: Periodic BC
Hi again, is there a way to implement periodic boundry conditions in fipy? for instance: field(x=L)==field(x=-L) Thanks Francisco El mar, 17-05-2016 a las 09:44 -0400, Daniel Wheeler escribió: > 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 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 ] > > > ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]