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