Re: Moving boundary in FiPy

[Zhekai Deng]

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?

Here are my boundary condition (robin)
dc/dz = a*local_shear_rate*(1-c)*c

where a is a constant, local_shear_rate is a FaceVariable but only has
value on z direction.

Previously, when working in a fixed domain, I was using statement like this:
phi.faceGrad.constrain([a*local_shear_rate*(1-phi.faceValue)*phi.faceValue],
where = mesh.facesTop | mesh.physicalFaces["interface"])

where the facesTop and "interface" are marked either by Fipy or Gmsh.

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)))
wrote:

> See these for example,
>
>  -
> https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-internal-boundary-conditions
>  -
> https://stackoverflow.com/questions/33109744/fipy-interior-conducting-boundary-conditions
>  - https://www.mail-archive.com/fipy@nist.gov/msg03674.html
>  - https://gist.github.com/guyer/a61d5adfa9a050eb970a
>
>
> On Fri, Jun 22, 2018 at 10:43 AM, Daniel Wheeler
>  wrote:
> > On Thu, Jun 21, 2018 at 3:39 PM, Zhekai Deng
> >  wrote:
> >> Hi Daniel,
> >>
> >> Thanks for the reply. I have a quick question regarding to varying
> >> coefficients approach. I have two boundary conditions needed to apply,
> one
> >> on the free surface and other one below the free surface with certain
> depth.
> >> Thus, these are two the boundaries that I am trying to move as the tank
> gets
> >> depleted. These are two straight line moving at a constant speed as the
> tank
> >> depletes. Previously, because I solve the equation using fixed domain,
> I use
> >> gmsh to draw and mark the boundary. Does this prevent me from varying
> >> coefficient approach?
> >
> > No, not at all. The boundary conditions can be reformulated as sources
> > that only apply in the region of the interface.
> >
> > --
> > 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 ]
>
___
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

[Zhekai Deng]

Hi Daniel,

Thanks for the reply. I have a quick question regarding to varying
coefficients approach. I have two boundary conditions needed to apply, one
on the free surface and other one below the free surface with certain
depth. Thus, these are two the boundaries that I am trying to move as the
tank gets depleted. These are two straight line moving at a constant speed
as the tank depletes. Previously, because I solve the equation using fixed
domain, I use gmsh to draw and mark the boundary. Does this prevent me from
varying coefficient approach?

I understand that re-draw mesh at every small interval approach could be
slow but I am more worry about maybe the error or legitimacy of my naïve
ideas.

Thanks for the help!

Best,

Zhekai

On Tue, Jun 19, 2018 at 2:56 PM Daniel Wheeler 
wrote:

> On Tue, Jun 19, 2018 at 11:49 AM, Zhekai Deng
>  wrote:
> >
> > My initial guess would be at every time step (or a small interval), I
> will
> > need to reconstruct my new mesh and interpolate the old solution to the
> new
> > mesh. Because the fluids is being depleted, number of cells in new mesh
> will
> > be reduced and consequently less computation time.  I wonder do you might
> > have a better suggestion?
>
> I don't. My view is that it would still be easier and more efficient
> to use an approach with a fixed mesh and vary the coefficients /
> properties based on the fluid position.
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Moving boundary in FiPy

[Zhekai Deng]

Hi Daniel,

Thank you very much for the reply! When I initially saw your response, I
thought maybe it would be hard to implement in FiPy. However, I think given
the easy usage of FiPy and my previous FiPy build, It might still be easier
for me to implement in FiPy. My geometry is pretty simple: a triangle
stitched to the bottom of a rectangular (cross section of a silo). I am
trying to solve an advection-diffusion equation in the fluid drainage in
a tank. The tank is initially filled with fluid which is depleted through
the bottom outlet when time starts. I know all the information within the
tank (velocity field, shape of free surface, and depletion rate), and I
also know how to solve my nonlinear advection diffusion equation in a fixed
boundary. My initial guess would be at every time step (or a small
interval), I will need to reconstruct my new mesh and interpolate the old
solution to the new mesh. Because the fluids is being depleted, number of
cells in new mesh will be reduced and consequently less computation time.
I wonder do you might have a better suggestion?

Best,

Zhekai

On Thu, May 10, 2018 at 1:13 PM Daniel Wheeler 
wrote:

> 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
>  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 ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Moving boundary in FiPy

[Zhekai Deng]

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 ]

Re: Memory leak using trilinos

[Zhekai Deng]

Hi Mike,

Thanks for your reply. I just tried the garbage collection and it seems did
not help very much with the leak issue for trilinos. My Pysparse throw me
Segmentation fault, which I don't really understand why.

For the check point approach, just to make sure that I understand it
correctly, are you using dump from Fipy.tools to store the solution
variable, and read them back ? I guess I also need to re-declare my
equation formulations, and boundary conditions each time I read them back.
Also, I wonder how to delete/remove all the previous object so that I could
have a "fresh" start after each check point iteration ?

Thanks,

Zhekai


On Tue, Jun 20, 2017 at 4:51 PM, Michael J. Waters 
wrote:

> Hi Zhekai Deng,
>
> I had a similar problem. In the end, I just created check point files and
> restarted frequently.
>
> I did find that garbage collecting helped slow the leak and,
> interestingly, also sped up code execution for me.
>
> Best,
>
> -Mike
>
> On 6/20/17 10:39 AM, Zhekai Deng wrote:
>
> HI All,
>
> I think I run into some memory leak issues using trilinos as my solver
> both for serial and parallel case. The total memory goes up from initial
> around 300 MB to the my system memory limit (like 29 GB), and program
> stops. I have attached my example code below.
>
> Somethings I am not sure if I implement it correctly:
> 1. in the example code, I basically solve multiple equations at the same
> time. I use list, and append my CellVarible together
>
> 2. Similar, I use list and append my equation together
>
> 3. Create coupled equation by looping through the list element, and solve
> this coupled equation.
>
> 4. During the solution step, I loop through the list element, and do
> updateOld()
>
> The results of the numerical calculation looks correct to me, I just have
> this memory issue that prevent me from going longer time steps because I
> run out of  memory.
>
> I wonder is there any workaround/fix to this problem ?
>
> I am using the build from Anaconda repository, default trillions solver,
> and my OS is Ubuntu 16.04.
>
>
> ___
> fipy mailing listfipy@nist.govhttp://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 ]

Memory leak using trilinos

[Zhekai Deng]

HI All,

I think I run into some memory leak issues using trilinos as my solver both
for serial and parallel case. The total memory goes up from initial around
300 MB to the my system memory limit (like 29 GB), and program stops. I
have attached my example code below.

Somethings I am not sure if I implement it correctly:
1. in the example code, I basically solve multiple equations at the same
time. I use list, and append my CellVarible together

2. Similar, I use list and append my equation together

3. Create coupled equation by looping through the list element, and solve
this coupled equation.

4. During the solution step, I loop through the list element, and do
updateOld()

The results of the numerical calculation looks correct to me, I just have
this memory issue that prevent me from going longer time steps because I
run out of  memory.

I wonder is there any workaround/fix to this problem ?

I am using the build from Anaconda repository, default trillions solver,
and my OS is Ubuntu 16.04.
from fipy import *
import numpy as np
import  matplotlib as matplotlib_plot
import scipy.interpolate
import resource

def make2dList(rows, cols):
a=[]
for row in xrange(rows): a += [[0]*cols]
return a



mesh = Gmsh2D('''
lc = 1e-2;
Point(1) = {-1, 0, 0, lc};
Point(2) = {0, 0, 0, lc};
Point(3) = {1, 0, 0, lc};
Point(4) = {0,-0.1,0,lc};
Ellipse(5) = {1,2,2,4};
Ellipse(6) = {4,2,2,3};
Line(7) = {3,2};
Line(8) = {2,1};
Circle(9) = {1,2,3};
Line Loop(9) = {6,7, 8,5};
Line Loop(10) = {-6,-5,9};
Plane Surface (11) = {9};
Plane Surface (12) = {10};
Physical Surface("top") = {11};
Physical Surface("bottom") = {12};
Physical Line("interface") = {5, 6};
  ''' % locals())
  
Np = 3
phi =[];
for species_i in range(Np-1):
species = CellVariable(mesh=mesh, hasOld = True, value=1.0/Np)
phi.append(species)

xFace,yFace=mesh.faceCenters
Epsilon = 0.1
velocityX_Top = (1-Epsilon*Epsilon)/(Epsilon*Epsilon)*(yFace + Epsilon*np.sqrt(1-xFace*xFace))
velocityY_Top = (1-Epsilon*Epsilon)/Epsilon*(xFace*yFace)/(np.sqrt(1-xFace*xFace))

velocityX_Bottom = yFace
velocityY_Bottom = -xFace


velocityVector = FaceVariable(mesh=mesh, rank=1)
topFaces = numerix.unique(mesh.cellFaceIDs[..., mesh.physicalCells["top"].value].flatten().filled())
bottomFaces = numerix.unique(mesh.cellFaceIDs[..., mesh.physicalCells["bottom"].value].flatten().filled())

velocityVector[0, topFaces] = velocityX_Top[..., topFaces]
velocityVector[1, topFaces] = velocityY_Top[..., topFaces]

velocityVector[0, bottomFaces] = velocityX_Bottom[..., bottomFaces]
velocityVector[1, bottomFaces] = velocityY_Bottom[..., bottomFaces]

local_shear_rate = FaceVariable(mesh=mesh, rank=1)
local_shear_rate[0, topFaces] = 0.
local_shear_rate[1, topFaces] = 1.0
local_shear_rate[0, bottomFaces] = 0.
local_shear_rate[1, bottomFaces] = 0.

# Velocity Field Plot
#viewer_velocity = MatplotlibVectorViewer(vars=velocityVector,title="velocity_vector")
#viewer_velocity.plot()

Pe = 5.

D = FaceVariable(mesh=mesh, rank=1)
D[0,topFaces] = 0.
D[1,topFaces] = Epsilon*Epsilon/Pe
D[0,bottomFaces] = 0.
D[1,bottomFaces] = 0.

Lambda = make2dList(Np, Np)
radius = np.linspace(0.001,0.003,num=Np)
beta = 0.26
R = 0.075
c = (1-pow(Epsilon,2))/(R*pow(Epsilon,3))

for i in range(Np):
for j in range(Np):
Lambda[i][j] = c*beta*2.0*np.minimum(radius[i],radius[j])*np.log(radius[i]/radius[j])

coupledEqn = 0
eq = []
Segregation_Term = [0]*(Np-1)
Last_Species = 1.0

for species_i in range(Np-1):
print species_i
for species_j in range(Np -1):
Segregation_Term[species_i] = Segregation_Term[species_i] + Lambda[species_i][species_j]*phi[species_j]
Last_Species = Last_Species - phi[species_j]

Segregation_Term[species_i] = Segregation_Term[species_i] + Lambda[species_i][Np-1]*Last_Species

eq.append(TransientTerm(var=phi[species_i]) \
+ VanLeerConvectionTerm(coeff = velocityVector + Epsilon*local_shear_rate*Segregation_Term[species_i].faceValue, var=phi[species_i])\
== DiffusionTerm(coeff=D,var=phi[species_i]))

phi[species_i].faceGrad.constrain([1.0/Epsilon*local_shear_rate*Pe*(phi[species_i]*Segregation_Term[species_i]).faceValue], where = mesh.facesTop | mesh.physicalFaces["interface"]) #boundary conditions at the top and bottom
Last_Species = 1.0

for species_i in range(Np-1):
coupledEqn = coupledEqn & eq[species_i]

timeStepDuration =  1/2000.
steps = 1100

k = 0
for step in range(steps):
res = 1000
while res > 1e-6:
res = coupledEqn.sweep(dt=timeStepDuration)
	print "MaxRSS:", resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
print res

for species_i in range(Np-1):
phi[species_i].updateOld()

print step
__

Re: Assign FaceVariable value based on submesh ID

[Zhekai Deng]

Hi Jon,

Thank you. It works perfectly.

Best,

Zhekai

On Sun, May 21, 2017 at 12:13 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> Zhekai -
>
> Thank you for the explanation. You're right, I answered a much simpler
> question than you were asking.
>
> The following should do what you want:
>
> velocityVector = fp.FaceVariable(mesh=mesh, rank=1)
> topFaces = fp.numerix.unique(mesh.cellFaceIDs[...,
> mesh.physicalCells["top"].value].flatten().filled())
> velocityVector[0, topFaces] = velocityX_Top[..., topFaces]
> velocityVector[1, topFaces] = velocityY_Top[..., topFaces]
>
> This obtains all of the faces that surround each of the cells in "top",
> rearranges them into a long list, removes any empty values (different cells
> can have different numbers of faces), and then removes duplicates. The
> result is a list of the faceIDs for all the faces that bound all of the
> cells in "top".
>
>
> Alternatively,
>
> topFaces = (mesh.physicalCells["top"].faceValue != 0.).value
>
> This obtains the face-evaluated average of whether a cell is in "top".
> This relies on the fact that True and False can be treated pretty
> interchangeably with 1. and 0.. The result is a boolean mask which is True
> for all the faces that bound all of the cells in "top".
>
> - Jon
>
> > On May 12, 2017, at 12:13 PM, Zhekai Deng  northwestern.edu> wrote:
> >
> > Hi Jon,
> >
> > Thanks for answering that. However, I am afraid I am not fully
> understand how "physicalFaces" field for gmsh mesh could help me resolve
> this one.
> >
> > The "physicalFaces" field, If I understand correctly, is giving me the
> Physical Line (not Physical Surface) that I specified in gmsh. Thus, this
> is giving me the faces index on that specific line I specified in gmsh.
> What I want to achieve is that is there any function that takes sub-domain
> name as input (in my case "Physical Surface("top")) give me all the faces
> in this domain (not just faces on one line).
> >
> > If I just use mesh.physicalFaces["Top"], it will just return error
> because "Top" in my case is a "Physical Surface" not a "Physical Line".
> >
> > Any help will be appreciated. Thanks.
> >
> > Zhekai
> >
> > On Fri, May 12, 2017 at 8:54 AM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
> > There is also a .physicalFaces field defined for a gmsh mesh.
> >
> > > On May 11, 2017, at 6:39 PM, Zhekai Deng  northwestern.edu> wrote:
> > >
> > > Hi All,
> > >
> > > I wonder is there any way to specify the FaceVariable based on the
> naming from gmsh ?  For example, I name my sub domain in gmsh as "top" and
> "bottom". In Fipy, I would like to do something similar to following:
> > >
> > > velocityX_Top = (1-Epsilon*Epsilon)/(Epsilon*Epsilon)*(yFace +
> Epsilon*np.sqrt(1-xFace*xFace))
> > > velocityY_Top = (1-Epsilon*Epsilon)/Epsilon*(
> xFace*yFace)/(np.sqrt(1-xFace*xFace))
> > > velocityX_Bottom = yFace
> > > velocityY_Bottom = -xFace
> > >
> > > velocityVector = FaceVariable(mesh=mesh.physicalCells["top"], rank=1)
> > > velocityVector[0] = velocityX_Top
> > > velocityVector[1] = velocityY_Top
> > >
> > > velocityVector = FaceVariable(mesh=mesh.physicalCells["bottom"],
> rank=1)
> > > velocityVector[0] = velocityX_Bottom
> > > velocityVector[1] = velocityY_Bottom
> > >
> > > I think it gives me error because the mesh.physicalCells is cell
> index, not face index. But is there any way I could do get all face index
> using sub domain ID?
> > >
> > > I have attached the complete and minimal code to demonstrate what I
> try to achieve.
> > >
> > > Thanks!
> > >
> > > Zhekai
> > >
> > > ___
> > > 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 ]

Re: Assign FaceVariable value based on submesh ID

[Zhekai Deng]

Hi Jon,

Thanks for answering that. However, I am afraid I am not fully understand
how "physicalFaces" field for gmsh mesh could help me resolve this one.

The "physicalFaces" field, If I understand correctly, is giving me the
Physical Line (not Physical Surface) that I specified in gmsh. Thus, this
is giving me the faces index on that specific line I specified in gmsh.
What I want to achieve is that is there any function that takes sub-domain
name as input (in my case "Physical Surface("top")) give me all the faces
in this domain (not just faces on one line).

If I just use mesh.physicalFaces["Top"], it will just return error because
"Top" in my case is a "Physical Surface" not a "Physical Line".

Any help will be appreciated. Thanks.

Zhekai

On Fri, May 12, 2017 at 8:54 AM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> There is also a .physicalFaces field defined for a gmsh mesh.
>
> > On May 11, 2017, at 6:39 PM, Zhekai Deng  northwestern.edu> wrote:
> >
> > Hi All,
> >
> > I wonder is there any way to specify the FaceVariable based on the
> naming from gmsh ?  For example, I name my sub domain in gmsh as "top" and
> "bottom". In Fipy, I would like to do something similar to following:
> >
> > velocityX_Top = (1-Epsilon*Epsilon)/(Epsilon*Epsilon)*(yFace +
> Epsilon*np.sqrt(1-xFace*xFace))
> > velocityY_Top = (1-Epsilon*Epsilon)/Epsilon*(
> xFace*yFace)/(np.sqrt(1-xFace*xFace))
> > velocityX_Bottom = yFace
> > velocityY_Bottom = -xFace
> >
> > velocityVector = FaceVariable(mesh=mesh.physicalCells["top"], rank=1)
> > velocityVector[0] = velocityX_Top
> > velocityVector[1] = velocityY_Top
> >
> > velocityVector = FaceVariable(mesh=mesh.physicalCells["bottom"], rank=1)
> > velocityVector[0] = velocityX_Bottom
> > velocityVector[1] = velocityY_Bottom
> >
> > I think it gives me error because the mesh.physicalCells is cell index,
> not face index. But is there any way I could do get all face index using
> sub domain ID?
> >
> > I have attached the complete and minimal code to demonstrate what I try
> to achieve.
> >
> > Thanks!
> >
> > Zhekai
> >
> > ___
> > 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 ]

Assign FaceVariable value based on submesh ID

[Zhekai Deng]

Hi All,

I wonder is there any way to specify the FaceVariable based on the naming
from gmsh ?  For example, I name my sub domain in gmsh as "top" and
"bottom". In Fipy, I would like to do something similar to following:

velocityX_Top = (1-Epsilon*Epsilon)/(Epsilon*Epsilon)*(yFace +
Epsilon*np.sqrt(1-xFace*xFace))

velocityY_Top =
(1-Epsilon*Epsilon)/Epsilon*(xFace*yFace)/(np.sqrt(1-xFace*xFace))

velocityX_Bottom = yFace

velocityY_Bottom = -xFace


velocityVector = FaceVariable(mesh=mesh.physicalCells["top"], rank=1)

velocityVector[0] = velocityX_Top

velocityVector[1] = velocityY_Top


velocityVector = FaceVariable(mesh=mesh.physicalCells["bottom"], rank=1)

velocityVector[0] = velocityX_Bottom

velocityVector[1] = velocityY_Bottom


I think it gives me error because the mesh.physicalCells is cell index, not
face index. But is there any way I could do get all face index using sub
domain ID?


I have attached the complete and minimal code to demonstrate what I try to
achieve.


Thanks!


Zhekai
from fipy import *
import numpy as np
import  matplotlib as matplotlib_plot
import scipy.interpolate

mesh = Gmsh2D('''
lc = 2e-2;
Point(1) = {-1, 0, 0, lc};
Point(2) = {0, 0, 0, lc};
Point(3) = {1, 0, 0, lc};
Point(4) = {0,-0.1,0,lc};
Ellipse(5) = {1,2,2,4};
Ellipse(6) = {4,2,2,3};
Line(7) = {3,2};
Line(8) = {2,1};
Circle(9) = {1,2,3};
Line Loop(9) = {6,7, 8,5};
Line Loop(10) = {-6,-5,9};
Plane Surface (11) = {9};
Plane Surface (12) = {10};
Physical Surface("top") = {11};
Physical Surface("bottom") = {12};
Physical Line("interface") = {5, 6};
  ''' % locals())
  
phi_1 = CellVariable(name = "solution variable",
   mesh = mesh,
   hasOld = True,
   value = 0.5) # doctest: +GMSH

viewer = Matplotlib2DViewer(vars=phi_1, title="final solution", cmap =  
matplotlib_plot.cm.hot,datamin = 0., datamax = 1.)
#viewer.plot()
#
xFace,yFace=mesh.faceCenters
Epsilon = 0.1
velocityX_Top = (1-Epsilon*Epsilon)/(Epsilon*Epsilon)*(yFace + 
Epsilon*np.sqrt(1-xFace*xFace))
velocityY_Top = 
(1-Epsilon*Epsilon)/Epsilon*(xFace*yFace)/(np.sqrt(1-xFace*xFace))
velocityX_Bottom = yFace
velocityY_Bottom = -xFace
velocityVector = FaceVariable(mesh=mesh.physicalCells["top"], rank=1)
velocityVector[0] = velocityX_Top
velocityVector[1] = velocityY_Top___
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

[Zhekai Deng]

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 ?


Outlet_BC_value = Variable(value = 0.5, array = numerix.zeros((1,102)))

phi.constrain(Outlet_BC_value, where = mesh_2.exteriorFaces &
mesh_2.physicalFaces["Right"])


if for step in range(steps):

 Outlet_BC_value_new = scipy.interpolate.griddata(points, values,
(xFace_2[a], yFace_2[a]), method='nearest')

 Outlet_BC_value.setValue(Outlet_BC_value_new)

 #Then solve equation


Best,


Zhekai

On Tue, Oct 25, 2016 at 11:44 AM, Daniel Wheeler 
wrote:

> On Tue, Oct 25, 2016 at 12:34 PM, Zhekai Deng
>  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 ]
>
___
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

[Zhekai Deng]

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.

Best,

Zhekai

On Mon, Oct 24, 2016 at 9:00 AM, Daniel Wheeler 
wrote:

> On Mon, Oct 24, 2016 at 12:20 AM, Zhekai Deng
>  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 ]
>
___
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

[Zhekai Deng]

Sorry, I accidentally hit the "send" button. So I am just going to continue
with previous email. The problem with above approach is that while it
transfer the data successfully, the data is not in the right spatial
position. This could be because the order which fipy name in mesh_1 is
different from mesh_2. So this inlet_BC_value_new might have the position
index from mesh_1, but mesh_2 does not use the same index to order the
cell. Consequently, the data, when it is transferred to mesh_2, is at
different position.

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 ?

Also, the reason I initialized this Variable has 102 element is that I
checked on the boundary, there is 102 cells. I realized I f I don't
initialized with correct dimension, then I can't acquire data in the
solving step.

inlet_BC_value = Variable(value = 0.5, array = numerix.zeros((1,102)))


PS. Thanks Ian for providing the example before.


Any suggestion is appreciated.


Best,


Zhekai


On Sun, Oct 23, 2016 at 11:10 PM, Zhekai Deng <
zhekaideng2...@u.northwestern.edu> wrote:

> Hi All,
>
> I tried to implement this approach in my problem, and encountered some
> problems. To keep this post more general, I posted the pseudo code, and
> point out the problem. I could also send the actual code in private if
> someone is willing to take a look at it.
>
> Basically, I have two 2D gmsh created irregular mesh that connected to
> each other. I am trying to "copy" the value from the boundary of mesh_1,
> into the mesh_2. Please note that mesh_1 and mesh_2 both share this
> boundary. After playing around the gmsh, I found that I could draw two
> fields together. However, when I actually created the surface, I only
> created the surface corresponding to mesh_1 and mesh_2
>
>
> mesh_1 = Gmsh2D(''
>
> 
>
> Line Loop(9) = {6,7, 8,5};
>
> Line Loop(10) = {-6,-5,9};
>
> Plane Surface (11) = {9};
>
> Physical Surface("flowinglayer") = {11};
>
> Physical Line("Bottom_Left") = {5};
>
> Physical Line("Bottom_Right") = {6};
>
> ''' % locals())
>
>
> mesh_2 = Gmsh2D('''
>
> .
>
> Line Loop(9) = {6,7, 8,5};
>
> Line Loop(10) = {-6,-5,9};
>
> Plane Surface (12) = {10};
>
> Physical Surface("bulkflow") = {12};
>
> Physical Line("Bottom_Left") = {5};
>
> Physical Line("Bottom_Right") = {6};
>
> ''' % locals())
>
>
> As you could see, mesh_1 only corresponding to line loop 9, mesh_2 only
> corresponding to line loop 10. The reason I do this is it allows the mesh
> between the interface of the mesh_1 and mesh_2 to mesh properly. Then I
> created variable to transfer data, as suggested by Ian.
>
>
> inlet_BC_value = Variable(value = 0.5, array = numerix.zeros((1,102)))
>
> phi_2.constrain(inlet_BC_value, where = mesh_2.exteriorFaces &
> mesh_2.physicalFaces["Bottom_Right"])
>
>
>
> Finally, solve the equation and update this "inlet_BC_value"
>
>
> for step in range(steps):
>
> res_1 = 1000
>
> res_2 = 1000
>
> inlet_BC_value_new = phi_1.faceValue[(mesh_1.exteriorFaces &
> mesh_1.physicalFaces["Bottom_Right"]).value]
>
> inlet_BC_value.value = np.float64(inlet_BC_value_new)
>
>
> while np.maximum(res_1,res_2) > 1e-4:
>
> res_1 = eq_1.sweep(var=phi_1,dt=timeStepDuration)
>
> res_2 = eq_2.sweep(var=phi_2,dt=timeStepDuration) # doctest: +GMSH
>
>
> phi_1.updateOld()
>
> phi_2.updateOld()
>
> viewer_1.plot()
>
> viewer_2.plot()
>
> print step
>
>
>
>
> On Wed, Oct 12, 2016 at 11:14 PM, Zhekai Deng  northwestern.edu> wrote:
>
>> I see. Thanks for clarifying it ! I will implement this over the weekend
>> and to see how it works out. Thank you.
>>
>> Best,
>>
>> Zhekai
>>
>> On Wed, Oct 12, 2016 at 8:15 AM, Campbell, Ian <
>> i.campbel...@imperial.ac.uk> wrote:
>>
>>> Dear Zhekai,
>>>
>>>
>>>
>>> My October 11th post to this list affects my recommendation to you from
>>> October 4th. You’ll probably be fine if you aren’t running your simulation
>>> for long, but if you are running it for long periods, it will likely become
>>> cripplingly slow. I've implemented the corrected version below for you,
>>> which should not have a memory leak because of the method of updating the
>>> boundary condition.
>>>
>>>
>>>
>>> I hope that helps

Re: FW: how to set up data transfer between two adjacent nonuniform meshs

[Zhekai Deng]

Hi All,

I tried to implement this approach in my problem, and encountered some
problems. To keep this post more general, I posted the pseudo code, and
point out the problem. I could also send the actual code in private if
someone is willing to take a look at it.

Basically, I have two 2D gmsh created irregular mesh that connected to each
other. I am trying to "copy" the value from the boundary of mesh_1, into
the mesh_2. Please note that mesh_1 and mesh_2 both share this boundary.
After playing around the gmsh, I found that I could draw two fields
together. However, when I actually created the surface, I only created the
surface corresponding to mesh_1 and mesh_2


mesh_1 = Gmsh2D(''



Line Loop(9) = {6,7, 8,5};

Line Loop(10) = {-6,-5,9};

Plane Surface (11) = {9};

Physical Surface("flowinglayer") = {11};

Physical Line("Bottom_Left") = {5};

Physical Line("Bottom_Right") = {6};

''' % locals())


mesh_2 = Gmsh2D('''

.

Line Loop(9) = {6,7, 8,5};

Line Loop(10) = {-6,-5,9};

Plane Surface (12) = {10};

Physical Surface("bulkflow") = {12};

Physical Line("Bottom_Left") = {5};

Physical Line("Bottom_Right") = {6};

''' % locals())


As you could see, mesh_1 only corresponding to line loop 9, mesh_2 only
corresponding to line loop 10. The reason I do this is it allows the mesh
between the interface of the mesh_1 and mesh_2 to mesh properly. Then I
created variable to transfer data, as suggested by Ian.


inlet_BC_value = Variable(value = 0.5, array = numerix.zeros((1,102)))

phi_2.constrain(inlet_BC_value, where = mesh_2.exteriorFaces &
mesh_2.physicalFaces["Bottom_Right"])



Finally, solve the equation and update this "inlet_BC_value"


for step in range(steps):

res_1 = 1000

res_2 = 1000

inlet_BC_value_new = phi_1.faceValue[(mesh_1.exteriorFaces &
mesh_1.physicalFaces["Bottom_Right"]).value]

inlet_BC_value.value = np.float64(inlet_BC_value_new)


while np.maximum(res_1,res_2) > 1e-4:

res_1 = eq_1.sweep(var=phi_1,dt=timeStepDuration)

res_2 = eq_2.sweep(var=phi_2,dt=timeStepDuration) # doctest: +GMSH


phi_1.updateOld()

phi_2.updateOld()

viewer_1.plot()

viewer_2.plot()

print step




On Wed, Oct 12, 2016 at 11:14 PM, Zhekai Deng <
zhekaideng2...@u.northwestern.edu> wrote:

> I see. Thanks for clarifying it ! I will implement this over the weekend
> and to see how it works out. Thank you.
>
> Best,
>
> Zhekai
>
> On Wed, Oct 12, 2016 at 8:15 AM, Campbell, Ian <
> i.campbel...@imperial.ac.uk> wrote:
>
>> Dear Zhekai,
>>
>>
>>
>> My October 11th post to this list affects my recommendation to you from
>> October 4th. You’ll probably be fine if you aren’t running your simulation
>> for long, but if you are running it for long periods, it will likely become
>> cripplingly slow. I've implemented the corrected version below for you,
>> which should not have a memory leak because of the method of updating the
>> boundary condition.
>>
>>
>>
>> I hope that helps.
>>
>>
>>
>> Best regards,
>>
>>
>>
>> -  Ian
>>
>>
>>
>> phi_1_initial, phi_2_initial = 0.0, 0.0
>> D1, D2 = 1.5, 1.75
>>
>> dx_1 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
>> 0.16221167, 0.10838638, 0.03806023]
>> dx_2 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
>> 0.16221167, 0.10838638, 0.03806023]
>>
>> mesh_1 = Grid1D(dx=dx_1, Lx=1.0)
>> mesh_2 = Grid1D(dx=dx_2, Lx=1.0)
>>
>> phi_1 = CellVariable(mesh_1, value=phi_1_initial)
>> phi_2 = CellVariable(mesh_2, value=phi_2_initial)
>>
>> eq_1 = TransientTerm() == ExplicitDiffusionTerm(coeff=D1)
>> eq_2 = TransientTerm() == ExplicitDiffusionTerm(coeff=D2)
>>
>> phi_1.faceGrad.constrain(0.5, where=mesh_1.facesLeft)
>> phi_1.constrain(1.0, where=mesh_1.facesRight)
>> inlet_BC_value = Variable()
>> phi_2.faceGrad.constrain(inlet_BC_value, where=mesh_2.facesLeft)
>> phi_2.faceGrad.constrain(-0.5, where=mesh_2.facesRight)
>>
>> timeStepDuration = 0.9 * min(dx_1)**2 / (2 * max(D1, D2))
>> steps = 100
>> viewer = Viewer(vars=(phi_1, phi_2))
>>
>> for step in range(steps):
>> eq_1.solve(var=phi_1, dt = timeStepDuration)
>> inlet_BC_value_new = phi_1.faceValue[mesh_2.faceCenters == 1.0]
>> # Acquire the data
>> inlet_BC_value.setValue(float(inlet_BC_value_new))
>>   # Apply that data on the 2nd mesh
>> eq_2.solve(var=phi_2, dt = timeStepDuration)
>> viewer.plot()
>>
>>
>>
>>
>>
>> *From:* fipy-boun...@

Re: Solving convection only problem

[Zhekai Deng]

Thanks for the reply. I searched more and just realized this question has
been answered before. I then switched to VanLeerConvectionTerm, It turns
out this method is able to keep the shock pretty well which has been
mentioned in the previous posts. I might try CLAWPACK in the future, but I
think the VanLeer terms is good for my project for now. Thanks again.

Best,

Zhekai

On Thu, Oct 20, 2016 at 9:05 AM, Daniel Wheeler 
wrote:

> 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
>  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 ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Solving convection only problem

[Zhekai Deng]

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
from fipy import *
import numpy as np
import  matplotlib as matplotlib_plot

  
mesh_2 = Gmsh2D('''
lc = 1e-2;
Point(1) = {-1, 0, 0, lc};
Point(2) = {0, 0, 0, lc};
Point(3) = {1, 0, 0, lc};
Point(4) = {0,-0.1,0,lc};
Ellipse(5) = {1,2,2,4};
Ellipse(6) = {4,2,2,3};
Line(7) = {3,2};
Line(8) = {2,1};
Circle(9) = {1,2,3};
Line Loop(9) = {6,7, 8,5};
Line Loop(10) = {-6,-5,9};
Plane Surface (12) = {10};
Physical Surface("bulk_flow") = {12};
Physical Line("Bottom_Left") = {5};
Physical Line("Bottom_Right") = {6};
  ''' % locals())  
  

phi_2 = CellVariable(name = "solution variable",
   mesh = mesh_2,
   value = 0.) # doctest: +GMSH

xFace_2,yFace_2=mesh_2.faceCenters


velocityX_2 = yFace_2
velocityY_2 = -xFace_2
velocityVector_2 = FaceVariable(mesh=mesh_2, rank=1)
velocityVector_2[0] = velocityX_2
velocityVector_2[1] = velocityY_2

exit_velocityX_2 = yFace_2
exit_velocityY_2 = -xFace_2
exit_velocityVector_2 = FaceVariable(mesh=mesh_2, rank=1)
exit_velocityVector_2[0] = exit_velocityX_2
exit_velocityVector_2[1] = exit_velocityY_2

out_flow_mask_2 = ((xFace_2 <= 0.) & (yFace_2 < 0.))
exteriorCoeff_2 = FaceVariable(mesh_2, value=exit_velocityVector_2, rank=1)
exteriorCoeff_2.setValue([[0.],[0.]], where=  ~(mesh_2.exteriorFaces & 
mesh_2.physicalFaces["Bottom_Left"]))
velocityVector_2.setValue([[0.],[0.]], where=  mesh_2.exteriorFaces & 
mesh_2.physicalFaces["Bottom_Left"])

 
eq = TransientTerm() \
+ PowerLawConvectionTerm(coeff = velocityVector_2)\
+ ImplicitSourceTerm(exteriorCoeff_2.divergence) == 0

phi_2.constrain(1., where = mesh_2.exteriorFaces & 
mesh_2.physicalFaces["Bottom_Right"]) # doctest: +GMSH


viewer_2 = Matplotlib2DViewer(vars=phi_2, title="final solution", cmap =  
matplotlib_plot.cm.hot,datamin = 0., datamax = 1.)

timeStepDuration =  1/200.
steps = 800
for step in range(steps):
eq.sweep(var=phi_2,dt=timeStepDuration) # doctest: +GMSH
viewer_2.plot()___
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

[Zhekai Deng]

I see. Thanks for clarifying it ! I will implement this over the weekend
and to see how it works out. Thank you.

Best,

Zhekai

On Wed, Oct 12, 2016 at 8:15 AM, Campbell, Ian 
wrote:

> Dear Zhekai,
>
>
>
> My October 11th post to this list affects my recommendation to you from
> October 4th. You’ll probably be fine if you aren’t running your simulation
> for long, but if you are running it for long periods, it will likely become
> cripplingly slow. I've implemented the corrected version below for you,
> which should not have a memory leak because of the method of updating the
> boundary condition.
>
>
>
> I hope that helps.
>
>
>
> Best regards,
>
>
>
> -  Ian
>
>
>
> phi_1_initial, phi_2_initial = 0.0, 0.0
> D1, D2 = 1.5, 1.75
>
> dx_1 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
> 0.16221167, 0.10838638, 0.03806023]
> dx_2 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
> 0.16221167, 0.10838638, 0.03806023]
>
> mesh_1 = Grid1D(dx=dx_1, Lx=1.0)
> mesh_2 = Grid1D(dx=dx_2, Lx=1.0)
>
> phi_1 = CellVariable(mesh_1, value=phi_1_initial)
> phi_2 = CellVariable(mesh_2, value=phi_2_initial)
>
> eq_1 = TransientTerm() == ExplicitDiffusionTerm(coeff=D1)
> eq_2 = TransientTerm() == ExplicitDiffusionTerm(coeff=D2)
>
> phi_1.faceGrad.constrain(0.5, where=mesh_1.facesLeft)
> phi_1.constrain(1.0, where=mesh_1.facesRight)
> inlet_BC_value = Variable()
> phi_2.faceGrad.constrain(inlet_BC_value, where=mesh_2.facesLeft)
> phi_2.faceGrad.constrain(-0.5, where=mesh_2.facesRight)
>
> timeStepDuration = 0.9 * min(dx_1)**2 / (2 * max(D1, D2))
> steps = 100
> viewer = Viewer(vars=(phi_1, phi_2))
>
> for step in range(steps):
> eq_1.solve(var=phi_1, dt = timeStepDuration)
> inlet_BC_value_new = phi_1.faceValue[mesh_2.faceCenters == 1.0] #
> Acquire the data
> inlet_BC_value.setValue(float(inlet_BC_value_new))
>   # Apply that data on the 2nd mesh
> eq_2.solve(var=phi_2, dt = timeStepDuration)
> viewer.plot()
>
>
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Zhekai Deng
> *Sent:* 05 October 2016 05:20
> *To:* fipy@nist.gov
> *Subject:* Re: how to set up data transfer between two adjacent
> nonuniform meshs
>
>
>
> Hi Ian,
>
>
>
> Thank you very much for the example and pointing out the additional
> reference. Those are very useful informations. This approach seems
> intuitive and easy to incorporate into the code. I will try to incorporate
> it into my own code.  Thanks again.
>
>
>
> Best,
>
>
>
> Zhekai
>
>
>
> On Tue, Oct 4, 2016 at 5:23 AM, Campbell, Ian 
> wrote:
>
> Hi Zhekai,
>
>
>
> There is.
>
>
>
> Try the following, where as you suggest, the interpolated value at the
> right-most face of mesh_1 is used. That (outward flowing) flux value is
> then copied to the inlet boundary condition (Neumann/flux) for mesh_2. A
> diffusion equation with a different coefficient for your solution variable
> phi is defined on the 2nd mesh – this should work if you change the
> equation further. Both meshes are non-uniform and of unit length.
>
>
>
> I don’t know from your question what type or magnitude boundary conditions
> you’re interested in on the other faces, but nonetheless, you can see in
> the figures produced from the code below that the outlet flux on the
> right-most face of mesh_1 is equal to the inlet flux on the left-most face
> of mesh_2. There’s likely a cleaner way to do this, by updating the value
> of the boundary condition for mesh_2 rather than re-defining it, but this
> is a start.
>
>
>
> An additional reference for you, here is Dr. Guyer’s suitably-named Gist:
> https://gist.github.com/guyer/bb199559c00f6047d466daa18554d83d
>
>
>
> Let us know if that works for you.
>
>
>
> Best,
>
>
>
> -  Ian
>
>
>
> phi_1_initial, phi_2_initial = 0.0, 0.0
>
> D1, D2 = 1.5, 1.75
>
> inlet_BC_value_init = 2.0
>
>
>
> dx_1 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
> 0.16221167, 0.10838638, 0.03806023]
>
> dx_2 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
> 0.16221167, 0.10838638, 0.03806023]
>
>
>
> mesh_1 = Grid1D(dx=dx_1, Lx=1.0)
>
> mesh_2 = Grid1D(dx=dx_2, Lx=1.0)
>
>
>
> phi_1 = CellVariable(mesh_1, value=phi_1_initial)
>
> phi_2 = CellVariable(mesh_2, value=phi_2_initial)
>
>
>
> eq_1 = TransientTerm() == ExplicitDiffusionTerm(coeff=D1)
>
> eq_2 = TransientTerm() == ExplicitDiffusionTerm(coeff=D2)
>
>
>
> # Boundary Conditions
>
>

Re: how to set up data transfer between two adjacent nonuniform meshs

[Zhekai Deng]

Hi Ian,

Thank you very much for the example and pointing out the additional
reference. Those are very useful informations. This approach seems
intuitive and easy to incorporate into the code. I will try to incorporate
it into my own code.  Thanks again.

Best,

Zhekai

On Tue, Oct 4, 2016 at 5:23 AM, Campbell, Ian 
wrote:

> Hi Zhekai,
>
>
>
> There is.
>
>
>
> Try the following, where as you suggest, the interpolated value at the
> right-most face of mesh_1 is used. That (outward flowing) flux value is
> then copied to the inlet boundary condition (Neumann/flux) for mesh_2. A
> diffusion equation with a different coefficient for your solution variable
> phi is defined on the 2nd mesh – this should work if you change the
> equation further. Both meshes are non-uniform and of unit length.
>
>
>
> I don’t know from your question what type or magnitude boundary conditions
> you’re interested in on the other faces, but nonetheless, you can see in
> the figures produced from the code below that the outlet flux on the
> right-most face of mesh_1 is equal to the inlet flux on the left-most face
> of mesh_2. There’s likely a cleaner way to do this, by updating the value
> of the boundary condition for mesh_2 rather than re-defining it, but this
> is a start.
>
>
>
> An additional reference for you, here is Dr. Guyer’s suitably-named Gist:
> https://gist.github.com/guyer/bb199559c00f6047d466daa18554d83d
>
>
>
> Let us know if that works for you.
>
>
>
> Best,
>
>
>
> -  Ian
>
>
>
> phi_1_initial, phi_2_initial = 0.0, 0.0
>
> D1, D2 = 1.5, 1.75
>
> inlet_BC_value_init = 2.0
>
>
>
> dx_1 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
> 0.16221167, 0.10838638, 0.03806023]
>
> dx_2 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
> 0.16221167, 0.10838638, 0.03806023]
>
>
>
> mesh_1 = Grid1D(dx=dx_1, Lx=1.0)
>
> mesh_2 = Grid1D(dx=dx_2, Lx=1.0)
>
>
>
> phi_1 = CellVariable(mesh_1, value=phi_1_initial)
>
> phi_2 = CellVariable(mesh_2, value=phi_2_initial)
>
>
>
> eq_1 = TransientTerm() == ExplicitDiffusionTerm(coeff=D1)
>
> eq_2 = TransientTerm() == ExplicitDiffusionTerm(coeff=D2)
>
>
>
> # Boundary Conditions
>
> phi_1.faceGrad.constrain(0.5, where=mesh_1.facesLeft)
>
> phi_1.constrain(1.0, where=mesh_1.facesRight)
>
> phi_2.faceGrad.constrain(inlet_BC_value_init, where=mesh_2.facesLeft)
>
> phi_2.faceGrad.constrain(-0.5, where=mesh_2.facesRight)
>
>
>
> timeStepDuration = 0.9 * min(dx_1)**2 / (2 * max(D1, D2))
>
> steps = 100
>
> viewer = Viewer(vars=(phi_1, phi_2))
>
>
>
> for step in range(steps):
>
> eq_1.solve(var=phi_1, dt = timeStepDuration)
>
> inlet_BC_value = phi_1.faceValue[mesh_2.faceCenters == 1.0] #
> Acquire the data
>
> phi_2.faceGrad.constrain(inlet_BC_value, where=mesh_2.facesLeft)#
> Apply that data on the 2nd mesh
>
> eq_2.solve(var=phi_2, dt = timeStepDuration)
>
> viewer.plot()
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Zhekai Deng
> *Sent:* 04 October 2016 03:29
> *To:* fipy@nist.gov
> *Subject:* how to set up data transfer between two adjacent nonuniform
> meshs
>
>
>
> Hello All,
>
>
>
> I wonder is there any way to allow the data share on the two nonuniform
> mesh's interface ? To explain this, let's say I have two meshes, mesh_1 and
> mesh_2. Different equations govern the mesh_1 and mesh_2, and solution
> variable on both meshes is phi.  I would like the outlet (on the right
> face) on mesh_1 served as inlet(on the left face) on mesh_2.
>
>
>
> I tired to concatenate two mesh. If they are uniform, this works fine.
> However, if two does not share the exact the same face (for example, one is
> uniform, and another is nonuniform), there seems to be problem. .
>
>
>
> Thus, I wonder is there any way to "interpolate" the phi.facevalue on
> mesh_1 outlet face, and apply this interpolation value into the inlet of
> mesh_2?
>
>
>
>
>
> Best,
>
>
>
> Zhekai
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

how to set up data transfer between two adjacent nonuniform meshs

[Zhekai Deng]

Hello All,

I wonder is there any way to allow the data share on the two nonuniform
mesh's interface ? To explain this, let's say I have two meshes, mesh_1 and
mesh_2. Different equations govern the mesh_1 and mesh_2, and solution
variable on both meshes is phi.  I would like the outlet (on the right
face) on mesh_1 served as inlet(on the left face) on mesh_2.

I tired to concatenate two mesh. If they are uniform, this works fine.
However, if two does not share the exact the same face (for example, one is
uniform, and another is nonuniform), there seems to be problem. .

Thus, I wonder is there any way to "interpolate" the phi.facevalue on
mesh_1 outlet face, and apply this interpolation value into the inlet of
mesh_2?


Best,

Zhekai
___
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

[Zhekai Deng]

Hi James,

Thanks for providing this demo to illustrate the problem. I don't have any
particular ideas exactly why the initial value of psi helps to reduce the
error and why transient problem  "advect" out. However, I have some
findings that may help on this.

Finding # 1: I noticed you have
tried ExponentialConvectionTerm, PowerLawConvectionTerm,
CentralDifferenceConvectionTerm,
all of these give exponential error. However, if you tried
UpwindConvectionTerm, this will give the right result on the steady state
solution. Thus, maybe using upwind method, the convection does not require
the value from downstream, consequently, the error from the downstream will
not excite the error toward the upstream. However, I am still very surprise
with the magnitude of the error from other methods, and how similar

Finding # 2: In the transient state problem, if I increase the mesh size
from 50*50 to 100*100. The error actually grows larger for the
ExponentialConvectionTerm, PowerLawConvectionTerm,
CentralDifferenceConvectionTerm.
To show this, if I change the mesh size to 100*100, the max psi value I
have is around 1.1 . However, if I change the mesh size to 50*50, the max
psi is 1.005366, which is several orders of magnitude lower in terms of
difference to the exact solution. This is also the case for the
UpwindConvectionTerm, however, the error for both mesh size are very small
(max(psi) = 1e-13 or 1e-14  + 1). So even in the transient state, the mesh
size appears to somehow amplify the error if we use finer mesh. I am
confused by this.

To now, it seems UpWindConvectionTerm appears to the the work around method
to this issue. Of course, I will be willing to learn more about what other
people think on this.

Best,

Zhekai




On Fri, Sep 16, 2016 at 8:53 AM, James Pringle  wrote:

> No worries -- I had to do it to figure out the problem in my more complex
> domain and equation... The issue which surprised me was that the value the
> variable was initialized to had an effect on the steady solution.
>
> Jamie
>
> On Fri, Sep 16, 2016 at 8:14 AM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
>
>> James -
>>
>> I think Daniel will have more insight into why this is happening and if
>> there is anything to be done about it besides artificial relaxation, but I
>> just want to say how much I appreciate your putting this together. This is
>> a very lucid illustration.
>>
>> - Jon
>>
>> > On Sep 15, 2016, at 5:13 PM, James Pringle  wrote:
>> >
>> > Dear FiPy users --
>> >
>> >This is a simple example of how and why fipy may fail to solve a
>> >steady advection diffusion problem, and how solving the transient
>> >problem can reduce the error. I also found something that was a
>> >surprise -- the "initial" condition of a steady problem can affect
>> >the solution for some solvers.
>> >
>> >At the end are two interesting questions for those who want to
>> >understand what FiPy is actually doing I admit to being a bit
>> >lost
>> >
>> >The equation I am solving is
>> >
>> > \Del\dot (D\del psi + u*psi) =0
>> >
>> >Where the diffusion D is 1, and u is a vector (1,0) -- so there is
>> >only a flow of speed -1 in the x direction.  I solve this equation
>> >on a 10 by 10 grid. The boundary conditions are no normal gradient
>> >on the y=0 and y=10 boundary:
>> >
>> > psi_y =0 at y=0 and y=10
>> >
>> >For the x boundary, I impose a value of x=1 on the inflow boundary
>> at x=10
>> >(this is a little tricky -- the way the equation is written, u is
>> >the negative of velocity).
>> >
>> > psi=1 at x=10
>> >
>> >and a no-normal-gradient condition at x=0.
>> >
>> > psi_x=0 at x=0
>> >
>> >since all of the domain and boundary is symmetrical with respect to
>> >y, we can assume psi_y=0 is zero everywhere. This reduces the
>> equation to
>> >
>> > psi_xx + psi_x =0
>> >
>> >The general solution to this equation is
>> >
>> > psi=C1+C2*exp(-x)
>> >
>> >Where C1 and C2 are constants. For these boundary conditions, C1=1
>> >and C2=0, so psi=1 everywhere.
>> >
>> >Now run the code SquareGrid_HomemadeDelaunay and look at figure(3)
>> >-- this is the plot of psi versus x, and you can see that it does
>> >not match the true solution of psi=1 everywhere! Instead, it
>> >appears to be decaying exponential. The blue line is actually just
>> >(1+exp(-x)). What is going on? It appears that small errors in the
>> >boundary condition at x=0 are allowing C2 to not be exactly 0, and
>> >this error is this exponential mode. The error is the artificial
>> >exiting of a correct mode of the interior equation, albeit one that
>> >should not be excited by these BC's.
>> >
>> >Interestingly, the size of this error depends on the value the psi
>> >is initially set to. If the line
>> >
>> >psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
>> >
>> >is cha

Re: Solving advection, segregation and diffusion equation

[Zhekai Deng]

Thanks for the reply. Those are very useful suggestions. I have some new
findings which may help shedding some light on this problem.

I have added TransientTerm and tried to see how the solution is developed
over time. I also remove the datamin and datamax input on the viewer, so
that I can see the full range of the solution variable data.

When I see the temporal development of the solution (attached in the
email), It looks stable. However, I noticed that the right outlet boundary
has extreme high concentration and did not outflow the flow domain. I think
It has due to the improper setup of my right outlet condition. Since the
materials keep building up on the right outlet boundary, it appears to
affect the solution in the actual flow domain(also attached in email) and
drives the solution unstable. This could be the main reason that solving
steady state equation directly never converges.

I then begin to look into how to properly set up "open outlet" boundary,
the closest I could find is following:

https://www.mail-archive.com/fipy@nist.gov/msg02797.html
and in the section of "Applying fixed flux boundary conditions" of official
website
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html

I tried to adapt the official method listed on the official website, but no
success is achieved.  Thus, I wonder is it possible to point out some
examples that used this method to define fixed flux outlet boundary
condition ? You may find my current progress on this part in the code I
attached.

Another thing I noticed is that,  there is a reported issue on robin
boundary condition:

https://github.com/usnistgov/fipy/issues/426

Would it be possible to explain more about this ? I don't think the up &
bottom boundary condition in my case causes the non convergent problem,
because from the Transient equation, the segregation of the material starts
to develop at the top&bottom boundary, which is consistent with my
expectation.

I have attached my current progress of the code, and the temporal
development of the solution.

Best,

Zhekai






On Thu, Sep 15, 2016 at 9:33 AM, Daniel Wheeler 
wrote:

> On Wed, Sep 14, 2016 at 7:52 PM, Guyer, Jonathan E. Dr. (Fed)
>  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 ]
>
from fipy import *

# mesh set up
nx = 100.
nz = 100.
dx = 1/nx
dz = 1/nz

mesh = Grid2D(dx=dx, dy=dz, nx=nx, ny=nz)
phi = CellVariable(name = "solution variable",mesh = mesh, value = 0.5)


# the follwing set up the velocity and shear rate in vector form.

# for those direction that I don't need it, I set its value as zero.

velocityX = mesh.faceCenters[1] # linear velocity profile
velocityY = FaceVariable(mesh=mesh, value = 0) # zero velocity in y direction
velocityVector = FaceVariable(mesh=mesh, rank=1)
velocityVector[0] = velocityX
velocityVector[1] = velocityY


local_shear_rate_X = FaceVariable(mesh=mesh, value = 0) # this term only acts 
on normal direction
local_shear_rate_Y = FaceVariable(mesh=mesh, value = 1) # linear velocity 
profile du/dz
local_shear_rate = FaceVariable(mesh=mesh, rank=1)
local_shear_rate[0] = local_shear_rate_X
local_shear_rate[1] = local_shear_rate_Y

Pe = 100. # some constants

#Setting up the diffusion part

D_Matrix = [[[0, 0],

 [0, 1/Pe]]]

#This is the new part to account for outflow fixed flux boundary (this does not 
set up correctly)
outflowflux = FaceVariable(mesh = mesh, value = phi.faceValue*velocityVector)
velocityVector.constrain(0.,mesh.facesRight)

# equation set up
eq = (TransientTerm() + UpwindConvectionTerm(coeff = velocityVector)\
+ UpwindConvectionTerm(coeff =local_shear_rate*(1-phi).faceValue)\
- (mesh.facesRight * outflowflux).divergence \
== DiffusionTerm (coeff=D_Matrix))


# set up the boundary conditions

phi.constrain(0.5, where = mesh.facesLeft) # fixed the inlet as 0.5 
concentration

phi.faceGrad.constrain(local_shear_rate*Pe*(1-phi.faceValue)*phi.faceValue, 
where = mesh.facesTop | mesh.facesBottom) #boundary conditions at the top and 
bottom


#phi.constrain(0, where = mesh.facesRight) #not working outlet boundary 
condition
#phi.faceGrad.constrain(velocityVector*phi.faceValue,where = mesh.facesRight) # 
problematic outlet boundary conditions

# Viewer set up

viewer = Matplotlib2

Re: Solving advection, segregation and diffusion equation

[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).

Solving advection, segregation and diffusion equation

[Zhekai Deng]

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" that we want to constrain ?
2. In my problem specifically,  I would like to implement convection inlet
flux and convection outlet flux. This seems is against the Fipy default no
flux condition.  Any suggestion to on how to implement it ?

Best,

Zhekai
from fipy import *
# mesh set up
nx = 300.
nz = 300.
dx = 1/nx
dz = 1/nz

mesh = Grid2D(dx=dx, dy=dz, nx=nx, ny=nz)
phi = CellVariable(name = "solution variable",mesh = mesh, value = 0.5)

# the follwing set up the velocity and shear rate in vector form.
# for those direction that I don't need it, I set its value as zero.
velocityX = mesh.faceCenters[1] # linear velocity profile
velocityY = FaceVariable(mesh=mesh, value = 0) # zero velocity in y direction
velocityVector = FaceVariable(mesh=mesh, rank=1)
velocityVector[0] = velocityX
velocityVector[1] = velocityY

shear_rate_X = FaceVariable(mesh=mesh, value = 0) # linear velocity profile
shear_rate_Y = FaceVariable(mesh=mesh, value = 1) # zero velocity in y direction
shear_rate = FaceVariable(mesh=mesh, rank=1)
shear_rate[0] = shear_rate_X
shear_rate[1] = shear_rate_Y


Pe = 100. # some constants

#Setting up the diffusion part
Diffusion_X = FaceVariable(mesh=mesh, value = 0)
Diffusion_Y = FaceVariable(mesh=mesh, value = 1/Pe)
D_Vector = FaceVariable(mesh=mesh, rank=1)
D_Vector[0] = Diffusion_X
D_Vector[1] = Diffusion_Y

# mask represents where the boundary condition is applied
mask = mesh.facesTop | mesh.facesBottom 
phi.constrain(0.5, where = mesh.facesLeft) # fixed the inlet as 0.5 
concentration
phi.faceGrad.constrain([shear_rate*Pe*(1-phi.faceValue)*phi.faceValue],mask) 
#boundary conditions at the top and bottom (Not sure whether it is correct or 
not)

# Viewer set up
viewer = Matplotlib2DViewer(vars=phi,datamin=0.,datamax = 1., title="final 
solution")

# Begin the calculation of the pde
res = 1e+10
while res >1e-6: #wait until it converges
eq = (ExponentialConvectionTerm(coeff = velocityVector) + 
ExponentialConvectionTerm(coeff = shear_rate*(1-phi).faceValue) == 
DiffusionTerm (coeff=D_Vector))
res = eq.sweep(var = phi)
viewer.plot()
print res
#I implmented both advection and segregation term using 
ExponentialConvectionTerm. 
#However, I would want there is a inlet flux on the left, and outlet flux on 
the right due to only advection.
#Any suggestions on how to implement is very appreciated! 


matlab_demo.m
Description: application/vnd.wolfram.mathematica.package


equation_explain.pdf
Description: Adobe PDF document
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: mesh size changes 2D advection diffusion problem solution

[Zhekai Deng]

Thank for you taking the time clarify this to me. I implemented your
suggestions, and it works right now. I understand the TransientTerm part
now. And I have  a quick follow up question regarding to velocity field
part: So in my original problem, the velocity field is linear in x
direction, and no velocity in y direction. So I think using
"mesh.faceCenters" gives me linear velocity on both direction. I think I
get the idea why we should use faceCenters instead of mesh.y or mesh.x, so
I made following changes to my definition of velocity field by declaring
them directly as FaceVariable. Is this the right things to do ?

Best,

Zhekai




velocityX = FaceVariable(mesh=mesh, value = mesh.faceCenters[1]) # linear
velocity profile

velocityY = FaceVariable(mesh=mesh, value = 0) # zero velocity in y
direction

velocityVector = FaceVariable(mesh=mesh, rank=1)

velocityVector[0] = velocityX

velocityVector[1] = velocityY





for step in range(steps):

eq = (TransientTerm() == DiffusionTerm(coeff=D)

- ExponentialConvectionTerm(coeff = velocityVector))

eq.solve(var = phi, dt = timeStepDuration)

viewer.plot()


-



On Wed, Mar 9, 2016 at 8:35 AM, Daniel Wheeler 
wrote:

> Hi Zhekai,
>
> Try adding a `TransientTerm()` to the equation and just using the
> `faceCenters` as the velocity so it looks like
>
> eq = (TransientTerm() == DiffusionTerm(coeff=D) -
> ExponentialConvectionTerm(coeff = mesh.faceCenters))
>
> and making the time step large. I used 1e4 to test.  Using the
> faceCenters ensures that the inlet/outlet condition has the same value
> as the mesh is resized. FiPy doesn't interpolate those values near the
> boundary from cell to face. I think that using a TransientTerm() makes
> the solution unique. It's probably always a good idea to use a
> TransientTerm() anyway and it never hurts to do so.
>
> Cheers,
>
> Daniel
>
>
> On Tue, Mar 8, 2016 at 1:04 PM, Zhekai Deng
>  wrote:
> > Dear all,
> >
> > I have a probably simple problem: I tried to solve the 2D
> > Convection-Diffusion problem with a velocity field that is spatially
> > dependent. I initialize my solution variable, phi, as 0.5. I found out
> that
> > If I don't implement the inlet boundary condition on the left (Dirichlet
> on
> > the left, on the line 24) and just impose velocity field, the solution
> > changes a lot when I increase mesh size. However, If I implement inlet
> > boundary condition, the solution seems stays the same as the mesh size
> > increases. I am not completely sure why mesh sizes will play a role in
> the
> > first case. Could someone clarify this problem to me ? Thank you very
> much.
> > I have copied and pasted the code in the following, and attached the
> source
> > code in the email.
> >
> > ---
> >
> > from fipy import *
> >
> > # when I change the nx from 400 to 200 and to 100 and ny from 400 to 200
> and
> > to 100.
> >
> > # The solutions changes a lot. I don't know exactly why the mesh would
> > change
> >
> > # the solution this much.
> >
> > nx = 200.
> >
> > ny = 200.
> >
> > dx = 1/nx
> >
> > dy = 1/ny
> >
> > L = 1.
> >
> > mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
> >
> > #initial phi is equal to 0.5 in all the space
> >
> > phi = CellVariable(name = "solution variable",mesh = mesh, value = 0.5)
> >
> >
> > # the follwing set up the velocity in vector form.
> >
> > # for those direction that I don't need it, I set its value as zero.
> >
> > velocityX = CellVariable(value = mesh.y,mesh=mesh, name=r'$u_x$') #
> linear
> > velocity profile
> >
> > velocityY = CellVariable(value = 0,mesh=mesh, name=r'$u_y$') # zero
> velocity
> > in y direction
> >
> > velocityVector = CellVariable(mesh=mesh, name=r'$\vec{u}$', rank=1)
> >
> > velocityVector[0] = velocityX
> >
> > velocityVector[1] = velocityY
> >
> > D = 0.1
> >
> >
> > #This is the inlet boundary conditions that I tried to implement.
> >
> > #phi.constrain(1, where = mesh.facesLeft)
> >
> > viewer = Matplotlib2DViewer(vars=phi, title="final solution")
> >
> >
> >
> > eq = (0 == DiffusionTerm(coeff=D)
> >
> > - ExponentialConvectionTerm(coeff = velocityVector.faceValue))
> >
> > eq.solve(var = phi)
> >
> > viewer.plot()
> >
> > 
> >
> >
> >
> > 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 ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

mesh size changes 2D advection diffusion problem solution

[Zhekai Deng]

Dear all,

I have a probably simple problem: I tried to solve the 2D
Convection-Diffusion problem with a velocity field that is spatially
dependent. I initialize my solution variable, phi, as 0.5. I found out that
If I don't implement the inlet boundary condition on the left (Dirichlet on
the left, on the line 24) and just impose velocity field, the solution
changes a lot when I increase mesh size. However, If I implement inlet
boundary condition, the solution seems stays the same as the mesh size
increases. I am not completely sure why mesh sizes will play a role in the
first case. Could someone clarify this problem to me ? Thank you very
much.  I have copied and pasted the code in the following, and attached the
source code in the email.

---

from fipy import *

# when I change the nx from 400 to 200 and to 100 and ny from 400 to 200
and to 100.

# The solutions changes a lot. I don't know exactly why the mesh would
change

# the solution this much.

nx = 200.

ny = 200.

dx = 1/nx

dy = 1/ny

L = 1.

mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

#initial phi is equal to 0.5 in all the space

phi = CellVariable(name = "solution variable",mesh = mesh, value = 0.5)


# the follwing set up the velocity in vector form.

# for those direction that I don't need it, I set its value as zero.

velocityX = CellVariable(value = mesh.y,mesh=mesh, name=r'$u_x$') # linear
velocity profile

velocityY = CellVariable(value = 0,mesh=mesh, name=r'$u_y$') # zero
velocity in y direction

velocityVector = CellVariable(mesh=mesh, name=r'$\vec{u}$', rank=1)

velocityVector[0] = velocityX

velocityVector[1] = velocityY

D = 0.1


#This is the inlet boundary conditions that I tried to implement.

#phi.constrain(1, where = mesh.facesLeft)

viewer = Matplotlib2DViewer(vars=phi, title="final solution")



eq = (0 == DiffusionTerm(coeff=D)

- ExponentialConvectionTerm(coeff = velocityVector.faceValue))

eq.solve(var = phi)

viewer.plot()




Thank you very much!

Best,

Zhekai
from fipy import *
# when I change the nx from 400 to 200 and to 100 and ny from 400 to 200 and to 
100. 
# The solutions changes a lot. I don't know exactly why the mesh would change
# the solution this much.
nx = 200. 
ny = 200. 
dx = 1/nx
dy = 1/ny
L = 1.
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
#initial phi is equal to 0.5 in all the space
phi = CellVariable(name = "solution variable",mesh = mesh, value = 0.5)

# the follwing set up the velocity in vector form.
# for those direction that I don't need it, I set its value as zero.
velocityX = CellVariable(value = mesh.y,mesh=mesh, name=r'$u_x$') # linear 
velocity profile
velocityY = CellVariable(value = 0,mesh=mesh, name=r'$u_y$') # zero velocity in 
y direction
velocityVector = CellVariable(mesh=mesh, name=r'$\vec{u}$', rank=1)
velocityVector[0] = velocityX
velocityVector[1] = velocityY
D = 0.1

#This is the inlet boundary conditions that I tried to implement. 
#phi.constrain(1, where = mesh.facesLeft)
viewer = Matplotlib2DViewer(vars=phi, title="final solution")


eq = (0 == DiffusionTerm(coeff=D) 
- ExponentialConvectionTerm(coeff = velocityVector.faceValue))
eq.solve(var = phi)
viewer.plot()___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]