Re: Solving advection, segregation and diffusion equation

2016-09-14 Thread Guyer, Jonathan E. Dr. (Fed)

> 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

2016-09-14 Thread Zhekai Deng
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

2016-09-14 Thread Daniel Wheeler
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

2016-09-14 Thread Daniel Wheeler
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

2016-09-14 Thread Francisco Vega Reyes

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 ]