Re: casting implicit Boundary Conditions in FiPy
On Thu, Jul 14, 2016 at 10:34 AM, Gopalakrishnan, Krishnakumar wrote: > Dear Dan, > > Thanks a lot for your reply. > > In the 'hacked' source term to handle this special boundary condition, i.e. > > fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * > mesh.facesRight).divergence)), > > my implicitCoeff turns out to be also a function of the x-coordinates of the > system. > > > > implicitCoeff = x * D * beta/((1 - beta*dx/2)) > > > > The fact that we are applying this source terms to the right face throws me > off, a bit. Usually, sources (i.e. generation terms) are defined at the > nodes, right ? Well, in the cell centred FV method everything is really "solved" at the cell centers or associated with the cell. The "implicitCoeff" is a face variable so it can change based on position. > > Consequently, the ‘x’ in my implicitCoeff is x = mesh.cellCenters[0], right > ? (And, NOT x = mesh.FaceCenters[0]) Looking back over the previous email, the "implicitCoeff" is neither a cell or face variable implicitCoeff = -D * k / (1. - k * dx / 2.) It's just a single value, I think. If you want to make that quantity depend on position then you can make it a face variable and also make it dependent on spatial positions of the faces. > The reason why I ask is, with your trick, we have zeroed out the diffusion > coefficient at the right face, and instead applying the trick you had > proposed by turning the boundary condition to the divergence of the > ImplicitSourceTerm applied to the right boundary by the Boolean > mesh.facesRight. Despite the fact that at first glance, ‘x’ looks like > mesh.FaceCenters[0], keeping with the concept of generation terms defined > only at node locations, the ‘x’ in my sourceterm coefficient must be x = > mesh.cellCenters[0]. The "implicitCoeff" needs to be defined at the faces as it's inside a divergence. It won't work otherwise. You may need to think carefully how to zero out the diffusion coefficients on the faces for 2D, but that can also be achieved based on spatial 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 ]
RE: casting implicit Boundary Conditions in FiPy
Dear Dan, Thanks a lot for your reply. In the 'hacked' source term to handle this special boundary condition, i.e. fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * mesh.facesRight).divergence)), my implicitCoeff turns out to be also a function of the x-coordinates of the system. implicitCoeff = x * D * beta/((1 - beta*dx/2)) The fact that we are applying this source terms to the right face throws me off, a bit. Usually, sources (i.e. generation terms) are defined at the nodes, right ? Consequently, the 'x' in my implicitCoeff is x = mesh.cellCenters[0], right ? (And, NOT x = mesh.FaceCenters[0]) The reason why I ask is, with your trick, we have zeroed out the diffusion coefficient at the right face, and instead applying the trick you had proposed by turning the boundary condition to the divergence of the ImplicitSourceTerm applied to the right boundary by the Boolean mesh.facesRight. Despite the fact that at first glance, 'x' looks like mesh.FaceCenters[0], keeping with the concept of generation terms defined only at node locations, the 'x' in my sourceterm coefficient must be x = mesh.cellCenters[0]. Is my understanding correct? Best Regards Krishna -Original Message- From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel WheelerTo: 'fipy@nist.gov' Subject: RE: casting implicit Boundary Conditions in FiPy I recommend trying both methods and implementing some 1D test cases to double check. I always do this no matter how much I trust the code. -Original Message- From: Gopalakrishnan, Krishnakumar Sent: Friday, July 1, 2016 7:44 PM To: 'fipy@nist.gov' Subject: RE: casting implicit Boundary Conditions in FiPy Hi Dan, In following up to the discussion we had so far, I have a question about a variant of the special-case boundary condition presented so far. The new (implicit Neumann) boundary condition for the 1D problem is now, d_phi/dx at boundary = alpha + beta*phi at boundary Where, alpha is just a constant. Using the trick you proposed, one can write as d_phi/dx at boundary = alpha + beta*(phi at last cell + (d_phi/dx at boundary * dx/2)) d_phi/dx at boundary = alpha/ (1 - beta*dx/2) + beta/((1 - beta*dx/2)) * phi _at_last cell That is, we have the boundary condition as the sum of an explicit and implicit term. It's very clear about what to do with the implicit term following your method and code. The implicit coefficient was used in the method fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * mesh.facesRight).divergence)). Now, my question is, for the explicit term in the BC, can't I just define in a regular Neumann BC style of definition ? , i.e. phi.facegrad.constrain([alpha/ (1 - beta*dx/2)], mesh.facesRight]) . i.e. will this ordinary Neumann correctly add up with the Implicit Boundary condition term , implemented using the ImplicitSourceTerm method ? Or, should I treat this explicit BC term as part of the PDE itself, (like the trick done for the implicit BC), ie. Do we need something along the lines of fp.SourceTerm((mesh.faceNormals * explicitCoeff* mesh.facesRight).divergence)) ? In either case, the diffusion coefficient along the problematic boundary has been forced to be zero. Regards Krishna implicitCoeff = beta/((1 - beta*dx/2)) explicitCoeff = lpha/ (1 - beta*dx/2) -Original Message- From: fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov> [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 23 June 2016 15:24 To: Multiple recipients of list mailto:fipy@nist.gov>> Subject: Re: casting implicit Boundary Conditions in FiPy On Fri, Jun 17, 2016 at 5:36 PM, Gopalakrishnan, Krishnakumar mailto:k.gopalakrishna...@imperial.ac.uk>> wrote: > > My problem models the solid diffusion in a spherical particle. Matter > diffuses from the centre of the particle and reacts at the surface. This > is captured in a normalised 1-D domain with suitable equations and > co-ordinate scaling. Particle-centre is represented at x = 0 > boundary, and surface of the sphere represented by x=1 boundary. > > > > Now, my meshing algorithm is a special one. My inter-nodal distance > (the thickness of each sub-shell of the sphere) is such that, all of > the internal shells have equal volume. This is done so that mass is > conserved within the domain. This practise is stemming from my finite > difference/finite element colleagues who advocate this. Hi Krishna, I'm not sure what the thinking is of your colleagues, but the size of the elements has little or no impact on conservation. In finite volume, the equations are conservative to at least the precision of the linear solver (if not more) independently from the accuracy of the solution. Of course there could be an issue with source terms that I
Re: casting implicit Boundary Conditions in FiPy
On Fri, Jul 1, 2016 at 2:44 PM, Gopalakrishnan, Krishnakumar wrote: > > Now, my question is, for the explicit term in the BC, can't I just define in > a regular Neumann BC style of definition ? , i.e. > phi.facegrad.constrain([alpha/ (1 - beta*dx/2)], mesh.facesRight]) . i.e. > will this ordinary Neumann correctly add up with the Implicit Boundary > condition term , implemented using the ImplicitSourceTerm method ? > > Or, should I treat this explicit BC term as part of the PDE itself, (like > the trick done for the implicit BC), ie. Do we need something along the lines > of fp.SourceTerm((mesh.faceNormals * explicitCoeff* > mesh.facesRight).divergence)) ? In either case, the diffusion > coefficient along the problematic boundary has been forced to be zero. I recommend trying both methods and implementing some 1D test cases to double check. I always do this no matter how much I trust the code. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
RE: casting implicit Boundary Conditions in FiPy
Hi Dan, In following up to the discussion we had so far, I have a question about a variant of the special-case boundary condition presented so far. The new (implicit Neumann) boundary condition for the 1D problem is now, d_phi/dx at boundary = alpha + beta*phi at boundary Where, alpha is just a constant. Using the trick you proposed, one can write as d_phi/dx at boundary = alpha + beta*(phi at last cell + (d_phi/dx at boundary * dx/2)) d_phi/dx at boundary = alpha/ (1 - beta*dx/2) + beta/((1 - beta*dx/2)) * phi _at_last cell That is, we have the boundary condition as the sum of an explicit and implicit term. It's very clear about what to do with the implicit term following your method and code. The implicit coefficient was used in the method fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * mesh.facesRight).divergence)). Now, my question is, for the explicit term in the BC, can't I just define in a regular Neumann BC style of definition ? , i.e. phi.facegrad.constrain([alpha/ (1 - beta*dx/2)], mesh.facesRight]) . i.e. will this ordinary Neumann correctly add up with the Implicit Boundary condition term , implemented using the ImplicitSourceTerm method ? Or, should I treat this explicit BC term as part of the PDE itself, (like the trick done for the implicit BC), ie. Do we need something along the lines of fp.SourceTerm((mesh.faceNormals * explicitCoeff* mesh.facesRight).divergence)) ? In either case, the diffusion coefficient along the problematic boundary has been forced to be zero. Regards Krishna implicitCoeff = beta/((1 - beta*dx/2)) explicitCoeff = lpha/ (1 - beta*dx/2) -Original Message- From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 23 June 2016 15:24 To: Multiple recipients of list Subject: Re: casting implicit Boundary Conditions in FiPy On Fri, Jun 17, 2016 at 5:36 PM, Gopalakrishnan, Krishnakumar wrote: > > My problem models the solid diffusion in a spherical particle. Matter > diffuses from the centre of the particle and reacts at the surface. This > is captured in a normalised 1-D domain with suitable equations and > co-ordinate scaling. Particle-centre is represented at x = 0 > boundary, and surface of the sphere represented by x=1 boundary. > > > > Now, my meshing algorithm is a special one. My inter-nodal distance > (the thickness of each sub-shell of the sphere) is such that, all of > the internal shells have equal volume. This is done so that mass is > conserved within the domain. This practise is stemming from my finite > difference/finite element colleagues who advocate this. Hi Krishna, I'm not sure what the thinking is of your colleagues, but the size of the elements has little or no impact on conservation. In finite volume, the equations are conservative to at least the precision of the linear solver (if not more) independently from the accuracy of the solution. Of course there could be an issue with source terms that I'm not seeing, but the diffusion and convection terms should be entirely conservative. > The problem is that, it is impossible to EXACTLY divide the shell into > integer number of iso-volume subshells. Thus, the scheme is chosen > such that, we get iso-volume shells for all inner sub-shells upto the last > shell. > The last shell's thickness (dx) is obviously and purposefully made > ultra-small following the discussions and Dan's solution proposed in this > thread. Clearly, this last shell's volume differs from the rest of the > inner subshells. > > > > Am I violating any conservation laws here by doing it ? I know that > finite volume is a conservative method. But this question nevertheless nags > me. I don't think so, what do the solutions show. Are the solutions conservative? > Is there any advantage to doing iso-volume subshells ( for the inner > shells) only to break this concept for the last shell ? I don't think so as long as the jump in the size of the cell volume is well within the precision of the solver. Three or four orders of magnitude should be okay. > Given that all the > solution dynamics happen at the surface (right boundary) , does a > simple geometric progression suffice with a very small dx at the right > boundary suffice ? Or is there any other optimal (structured) mesh > generation algorithm for choosing mesh sizes depending on > problem-type, that I need to refer to ? Don't know this one. Have you tried both methods and compared the solution? -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
> On Jun 23, 2016, at 10:23 AM, Daniel Wheeler > wrote: > > I'm not sure what the thinking is of your colleagues, but the size of > the elements has little or no impact on conservation. In finite > volume, the equations are conservative to at least the precision of > the linear solver (if not more) independently from the accuracy of the > solution. Of course there could be an issue with source terms that I'm > not seeing, but the diffusion and convection terms should be entirely > conservative. I think the bigger issue is that a Grid1D will not properly represent a spherical divergence, regardless of the volumes of the cells. I thought we had a set of SphericalGrid1D classes, but apparently not. It would be the same basic idea as CylindricalNonUniformGrid1D, but I think the cell volumes and face areas would need to be changed, I believe like this: def _calcFaceAreas(self): r = self._calcFaceCenters()[0] return r**2 def _calcCellVolumes(self): return super(CylindricalNonUniformGrid1D, self)._calcCellVolumes() / 3. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
On Fri, Jun 17, 2016 at 5:36 PM, Gopalakrishnan, Krishnakumar wrote: > > My problem models the solid diffusion in a spherical particle. Matter > diffuses from the centre of the particle and reacts at the surface. This > is captured in a normalised 1-D domain with suitable equations and > co-ordinate scaling. Particle-centre is represented at x = 0 boundary, and > surface of the sphere represented by x=1 boundary. > > > > Now, my meshing algorithm is a special one. My inter-nodal distance (the > thickness of each sub-shell of the sphere) is such that, all of the internal > shells have equal volume. This is done so that mass is conserved within the > domain. This practise is stemming from my finite difference/finite element > colleagues who advocate this. Hi Krishna, I'm not sure what the thinking is of your colleagues, but the size of the elements has little or no impact on conservation. In finite volume, the equations are conservative to at least the precision of the linear solver (if not more) independently from the accuracy of the solution. Of course there could be an issue with source terms that I'm not seeing, but the diffusion and convection terms should be entirely conservative. > The problem is that, it is impossible to EXACTLY divide the shell into > integer number of iso-volume subshells. Thus, the scheme is chosen such > that, we get iso-volume shells for all inner sub-shells upto the last shell. > The last shell's thickness (dx) is obviously and purposefully made > ultra-small following the discussions and Dan's solution proposed in this > thread. Clearly, this last shell's volume differs from the rest of the inner > subshells. > > > > Am I violating any conservation laws here by doing it ? I know that finite > volume is a conservative method. But this question nevertheless nags me. I don't think so, what do the solutions show. Are the solutions conservative? > Is there any advantage to doing iso-volume subshells ( for the inner shells) > only to break this concept for the last shell ? I don't think so as long as the jump in the size of the cell volume is well within the precision of the solver. Three or four orders of magnitude should be okay. > Given that all the > solution dynamics happen at the surface (right boundary) , does a simple > geometric progression suffice with a very small dx at the right boundary > suffice ? Or is there any other optimal (structured) mesh generation > algorithm for choosing mesh sizes depending on problem-type, that I need to > refer to ? Don't know this one. Have you tried both methods and compared the solution? -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
RE: casting implicit Boundary Conditions in FiPy
Hi Dan, Thanks for your reply, and for the tip to use Grid1D with variable dx. I was able to deploy it right away. There is a new complexity in this connection. My problem models the solid diffusion in a spherical particle. Matter diffuses from the centre of the particle and reacts at the surface. This is captured in a normalised 1-D domain with suitable equations and co-ordinate scaling. Particle-centre is represented at x = 0 boundary, and surface of the sphere represented by x=1 boundary. Now, my meshing algorithm is a special one. My inter-nodal distance (the thickness of each sub-shell of the sphere) is such that, all of the internal shells have equal volume. This is done so that mass is conserved within the domain. This practise is stemming from my finite difference/finite element colleagues who advocate this. The problem is that, it is impossible to EXACTLY divide the shell into integer number of iso-volume subshells. Thus, the scheme is chosen such that, we get iso-volume shells for all inner sub-shells upto the last shell. The last shell's thickness (dx) is obviously and purposefully made ultra-small following the discussions and Dan's solution proposed in this thread. Clearly, this last shell's volume differs from the rest of the inner subshells. Am I violating any conservation laws here by doing it ? I know that finite volume is a conservative method. But this question nevertheless nags me. Is there any advantage to doing iso-volume subshells ( for the inner shells) only to break this concept for the last shell ? Given that all the solution dynamics happen at the surface (right boundary) , does a simple geometric progression suffice with a very small dx at the right boundary suffice ? Or is there any other optimal (structured) mesh generation algorithm for choosing mesh sizes depending on problem-type, that I need to refer to ? Krishna -Original Message- From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 17 June 2016 14:50 To: Multiple recipients of list Subject: Re: casting implicit Boundary Conditions in FiPy On Thu, Jun 16, 2016 at 12:35 PM, Gopalakrishnan, Krishnakumar mailto:k.gopalakrishna...@imperial.ac.uk>> wrote: > Thanks. > > Yes, this Is indeed only first order accurate. I verified this by > successively cutting my dx by half, running your code, and comparing against > the Mathematica generated result. Each time dx is cut by half, the error is > also proportionately halved. > > This code requires me to take unreasonably small dx values. Since, neither > the solution, nor the gradients are available at each implicit time-steps, I > think that higher order schemes are probably ruled out. I think you're right. It's not easy. > I am thinking of using a variable mesh-sizing, let's say a log-spacing in 1D, > keeping a very fine spacing (ultra-small dx) for the last cell near the > boundary, and gradually taking bigger steps. This is also physically > consistent with my problem, wherein all the action takes place close to the > boundary, and nothing much is happening at the middle or left edges. Maybe it's possible to make the edge cell almost infinitely small and have the rest of the cells evenly spaced. > This brings me to another issue. I don't see a way to import a 1D .msh file > generated by gmsh into FiPy. I'm not sure if Gmsh does 1D meshes, you can always use a 2D mesh as a 1D mesh of course. > Secondly, I notice that there is an optimistic sounding grid1DBuilder > method. I couldn't find any help on how to use this method. Is this method > capable of creating the log-spaced mesh that I am considering ? Note that Grid1D takes an array for dx. So you can do >>> mesh = fp.Grid1D(dx=[0.5, 1.0, 2.0]) >>> print(mesh.cellCenters) [[ 0.25 1. 2.5 ]] As long as you can calculate the grid spacing then you don't need Gmsh. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov<mailto: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: casting implicit Boundary Conditions in FiPy
Dear Dr. Jonathan Guyer, Thank you very much for your help. For the benefit of the community wishing to use 1D meshes in gmsh format for any of their external applications, I am bringing into your attention a nifty script available from the 'fluidity' CFD project. For my 1D case, I have been able to use a very handy python script located here: https://github.com/FluidityProject/fluidity/blob/master/tools/interval.py The usage is as follows:. python interval.py --variable_dx=example.py 0.0 50.0 my_variable_mesh, will create a mesh in gmsh format called 'my_variabel_mesh.msh' in the current directory. The example.py can contain any arbitrary numpy function (with name val(X) ) for defining the mesh spacing. For example something like def val(X): return (5. / (X+1)) But clearly your script is much more amenable for a FiPy-centric workflow, since we don't have to worry about the intricacies of importing a 1D .msh file. Thank you for your help. Krishna -Original Message- From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Guyer, Jonathan E. Dr. (Fed) Sent: 17 June 2016 15:28 To: FIPY Subject: Re: casting implicit Boundary Conditions in FiPy Gmsh does do 1D meshes, and I've got experimental code that imports them, but it's not ready to merge, yet. In the meantime, this approach I've used for modeling semiconductor device contacts in 1D is probably better: n_thickness = 1e-6 # m p_thickness = 149e-6 # m grid_resolution = 5e-8 # m compression_factor = 0.8 compression_count = 10 compressed_dx = grid_resolution * compression_factor**numerix.arange(compression_count) compressed_length = compressed_dx.sum() compressed_dx = list(compressed_dx) n_uncompressed_thickness = n_thickness - 2 * compressed_length n_dx = n_uncompressed_thickness / round(n_uncompressed_thickness / grid_resolution) n_dx = compressed_dx[::-1] + [n_dx] * int(n_uncompressed_thickness / n_dx) + compressed_dx p_uncompressed_thickness = p_thickness - 2 * compressed_length p_dx = p_uncompressed_thickness / round(p_uncompressed_thickness / grid_resolution) p_dx = compressed_dx[::-1] + [p_dx] * int(p_uncompressed_thickness / p_dx + 0.5) + compressed_dx mesh = Grid1D(dx=n_dx + p_dx) > On Jun 17, 2016, at 9:50 AM, Daniel Wheeler > mailto:daniel.wheel...@gmail.com>> wrote: > > On Thu, Jun 16, 2016 at 12:35 PM, Gopalakrishnan, Krishnakumar > mailto:k.gopalakrishna...@imperial.ac.uk>> > wrote: >> Thanks. >> >> Yes, this Is indeed only first order accurate. I verified this by >> successively cutting my dx by half, running your code, and comparing against >> the Mathematica generated result. Each time dx is cut by half, the error is >> also proportionately halved. >> >> This code requires me to take unreasonably small dx values. Since, neither >> the solution, nor the gradients are available at each implicit time-steps, I >> think that higher order schemes are probably ruled out. > > I think you're right. It's not easy. > >> I am thinking of using a variable mesh-sizing, let's say a log-spacing in >> 1D, keeping a very fine spacing (ultra-small dx) for the last cell near the >> boundary, and gradually taking bigger steps. This is also physically >> consistent with my problem, wherein all the action takes place close to the >> boundary, and nothing much is happening at the middle or left edges. > > Maybe it's possible to make the edge cell almost infinitely small and > have the rest of the cells evenly spaced. > >> This brings me to another issue. I don't see a way to import a 1D .msh file >> generated by gmsh into FiPy. > > I'm not sure if Gmsh does 1D meshes, you can always use a 2D mesh as a > 1D mesh of course. > >> Secondly, I notice that there is an optimistic sounding grid1DBuilder >> method. I couldn't find any help on how to use this method. Is this method >> capable of creating the log-spaced mesh that I am considering ? > > Note that Grid1D takes an array for dx. So you can do > >>>> mesh = fp.Grid1D(dx=[0.5, 1.0, 2.0]) >>>> print(mesh.cellCenters) >[[ 0.25 1. 2.5 ]] > > As long as you can calculate the grid spacing then you don't need Gmsh. > > -- > Daniel Wheeler > > ___ > fipy mailing list > fipy@nist.gov<mailto: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<mailto: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: casting implicit Boundary Conditions in FiPy
Gmsh does do 1D meshes, and I've got experimental code that imports them, but it's not ready to merge, yet. In the meantime, this approach I've used for modeling semiconductor device contacts in 1D is probably better: n_thickness = 1e-6 # m p_thickness = 149e-6 # m grid_resolution = 5e-8 # m compression_factor = 0.8 compression_count = 10 compressed_dx = grid_resolution * compression_factor**numerix.arange(compression_count) compressed_length = compressed_dx.sum() compressed_dx = list(compressed_dx) n_uncompressed_thickness = n_thickness - 2 * compressed_length n_dx = n_uncompressed_thickness / round(n_uncompressed_thickness / grid_resolution) n_dx = compressed_dx[::-1] + [n_dx] * int(n_uncompressed_thickness / n_dx) + compressed_dx p_uncompressed_thickness = p_thickness - 2 * compressed_length p_dx = p_uncompressed_thickness / round(p_uncompressed_thickness / grid_resolution) p_dx = compressed_dx[::-1] + [p_dx] * int(p_uncompressed_thickness / p_dx + 0.5) + compressed_dx mesh = Grid1D(dx=n_dx + p_dx) > On Jun 17, 2016, at 9:50 AM, Daniel Wheeler wrote: > > On Thu, Jun 16, 2016 at 12:35 PM, Gopalakrishnan, Krishnakumar > wrote: >> Thanks. >> >> Yes, this Is indeed only first order accurate. I verified this by >> successively cutting my dx by half, running your code, and comparing against >> the Mathematica generated result. Each time dx is cut by half, the error is >> also proportionately halved. >> >> This code requires me to take unreasonably small dx values. Since, neither >> the solution, nor the gradients are available at each implicit time-steps, I >> think that higher order schemes are probably ruled out. > > I think you're right. It's not easy. > >> I am thinking of using a variable mesh-sizing, let's say a log-spacing in >> 1D, keeping a very fine spacing (ultra-small dx) for the last cell near the >> boundary, and gradually taking bigger steps. This is also physically >> consistent with my problem, wherein all the action takes place close to the >> boundary, and nothing much is happening at the middle or left edges. > > Maybe it's possible to make the edge cell almost infinitely small and > have the rest of the cells evenly spaced. > >> This brings me to another issue. I don't see a way to import a 1D .msh file >> generated by gmsh into FiPy. > > I'm not sure if Gmsh does 1D meshes, you can always use a 2D mesh as a > 1D mesh of course. > >> Secondly, I notice that there is an optimistic sounding grid1DBuilder >> method. I couldn't find any help on how to use this method. Is this method >> capable of creating the log-spaced mesh that I am considering ? > > Note that Grid1D takes an array for dx. So you can do > mesh = fp.Grid1D(dx=[0.5, 1.0, 2.0]) print(mesh.cellCenters) >[[ 0.25 1. 2.5 ]] > > As long as you can calculate the grid spacing then you don't need Gmsh. > > -- > 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: casting implicit Boundary Conditions in FiPy
On Thu, Jun 16, 2016 at 12:35 PM, Gopalakrishnan, Krishnakumar wrote: > Thanks. > > Yes, this Is indeed only first order accurate. I verified this by > successively cutting my dx by half, running your code, and comparing against > the Mathematica generated result. Each time dx is cut by half, the error is > also proportionately halved. > > This code requires me to take unreasonably small dx values. Since, neither > the solution, nor the gradients are available at each implicit time-steps, I > think that higher order schemes are probably ruled out. I think you're right. It's not easy. > I am thinking of using a variable mesh-sizing, let's say a log-spacing in 1D, > keeping a very fine spacing (ultra-small dx) for the last cell near the > boundary, and gradually taking bigger steps. This is also physically > consistent with my problem, wherein all the action takes place close to the > boundary, and nothing much is happening at the middle or left edges. Maybe it's possible to make the edge cell almost infinitely small and have the rest of the cells evenly spaced. > This brings me to another issue. I don't see a way to import a 1D .msh file > generated by gmsh into FiPy. I'm not sure if Gmsh does 1D meshes, you can always use a 2D mesh as a 1D mesh of course. > Secondly, I notice that there is an optimistic sounding grid1DBuilder > method. I couldn't find any help on how to use this method. Is this method > capable of creating the log-spaced mesh that I am considering ? Note that Grid1D takes an array for dx. So you can do >>> mesh = fp.Grid1D(dx=[0.5, 1.0, 2.0]) >>> print(mesh.cellCenters) [[ 0.25 1. 2.5 ]] As long as you can calculate the grid spacing then you don't need Gmsh. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
RE: casting implicit Boundary Conditions in FiPy
Thanks. Yes, this Is indeed only first order accurate. I verified this by successively cutting my dx by half, running your code, and comparing against the Mathematica generated result. Each time dx is cut by half, the error is also proportionately halved. This code requires me to take unreasonably small dx values. Since, neither the solution, nor the gradients are available at each implicit time-steps, I think that higher order schemes are probably ruled out. I am thinking of using a variable mesh-sizing, let's say a log-spacing in 1D, keeping a very fine spacing (ultra-small dx) for the last cell near the boundary, and gradually taking bigger steps. This is also physically consistent with my problem, wherein all the action takes place close to the boundary, and nothing much is happening at the middle or left edges. This brings me to another issue. I don't see a way to import a 1D .msh file generated by gmsh into FiPy. Secondly, I notice that there is an optimistic sounding grid1DBuilder method. I couldn't find any help on how to use this method. Is this method capable of creating the log-spaced mesh that I am considering ? Krishna -Original Message- From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 15 June 2016 17:50 To: Multiple recipients of list Subject: Re: casting implicit Boundary Conditions in FiPy On Wed, Jun 15, 2016 at 7:27 AM, Gopalakrishnan, Krishnakumar wrote: > > Dan. I was able validate that your code correctly implements the Implicit > Neumann Boundary Condition (in 1-D). > > Here is a link to the plot of the solution after a transient simulation for > 0.2 seconds. The plot on the left is generated in Mathematica, and that to > the right is that generated by the FiPy’s matplotlib viewer class. The two > are identical. > > https://imperialcollegelondon.box.com/s/lb8v1iqxb3rqhxsphn9cmhwcmh5jzp > iz Awesome that it worked. I think that you're right though regarding first versus second order boundary condition. I think it is first order in space. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
On Wed, Jun 15, 2016 at 7:27 AM, Gopalakrishnan, Krishnakumar wrote: > > Dan. I was able validate that your code correctly implements the Implicit > Neumann Boundary Condition (in 1-D). > > Here is a link to the plot of the solution after a transient simulation for > 0.2 seconds. The plot on the left is generated in Mathematica, and that to > the right is that generated by the FiPy’s matplotlib viewer class. The two > are identical. > > https://imperialcollegelondon.box.com/s/lb8v1iqxb3rqhxsphn9cmhwcmh5jzpiz Awesome that it worked. I think that you're right though regarding first versus second order boundary condition. I think it is first order in space. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
> On Jun 15, 2016, at 7:27 AM, Gopalakrishnan, Krishnakumar > wrote: > > Can we add this as a feature request in the project’s github page, perhaps > for a future release? By all means, please file issues on GitHub for any features you'd like or bugs you find. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
RE: casting implicit Boundary Conditions in FiPy
Hi, Dan. I was able validate that your code correctly implements the Implicit Neumann Boundary Condition (in 1-D). Here is a link to the plot of the solution after a transient simulation for 0.2 seconds. The plot on the left is generated in Mathematica, and that to the right is that generated by the FiPy’s matplotlib viewer class. The two are identical. https://imperialcollegelondon.box.com/s/lb8v1iqxb3rqhxsphn9cmhwcmh5jzpiz Can we add this as a feature request in the project’s github page, perhaps for a future release? Does Implicit Neumann BCs occur often in other disciplines ? Once again, thanks for your help in providing this solution approach. Krishna From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 10 June 2016 18:18 To: Multiple recipients of list Subject: Re: casting implicit Boundary Conditions in FiPy Hi Krishina and Ray, Thanks for the interesting discussion. I'm not 100% sure about everything that Krishina is asking for in the latter part of the discussion so I'm just going to address the code that Ray has developed below (my code is below). I think there is a way to handle the right hand side boundary condition both implicitly and with second order accuracy using an ImplicitSourceTerm boundary condition trick. Sort of similar to what is here, http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-fixed-flux-boundary-conditions. Like I said, I think the code below is both second order accurate (for fixed dx) and implicit. Extending this to 2D might raise a few issues. The grid spacing needs to be in the coefficient so, obviously, dx needs to be fixed in both x and y directions. Also, I'm not sure if I'm missing a factor of dx in the source term in 1D as I'm using both a divergence and an ImplicitSourceTerm so there is some question about volumes and face areas in 2D as well. I'm also confused about the signs, I had to flip the sign in front of the source to make it work. It seems to look right though in 1D and you can just take one massive time step to get to the answer. This will only work if k is negative otherwise it's unstable, right? The way I came up with the source was the following n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) and then solve for n.grad(phi) which gives n.grad(phi) = k * phi_P / (1 - dx * k / 2) where phi_P is the value of phi at the cell next to the boundary with the implicit boundary condition. Then we fake the outbound flux to be the expression on the right of the equation. Just as a general note it would be great in FiPy if we could come up with a nice way to write boundary conditions in scripts that did all these tricks implicitly without having to know all these background details about FV and how FiPy works. import fipy as fp nx = 50 dx = 1. mesh = fp.Grid1D(nx=nx, dx=dx) phi = fp.CellVariable(name="field variable", mesh=mesh, value=1.0) D = 1. k = -1. diffCoeff = fp.FaceVariable(mesh=mesh, value=D) diffCoeff.constrain(0., mesh.facesRight) valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain([k * phi.harmonicFaceValue], mesh.facesRight) # This is the problematic BC implicitCoeff = -D * k / (1. - k * dx / 2.) eq = (fp.TransientTerm() == fp.DiffusionTerm(diffCoeff) - \ fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * mesh.facesRight).divergence)) timeStep = 0.9 * dx**2 / (2 * D) timeStep = 10.0 steps = 800 viewer = fp.Viewer(vars=phi, datamax=1., datamin=0.) for step in range(steps): eq.solve(var=phi, dt=timeStep) viewer.plot() On Thu, Jun 9, 2016 at 12:02 PM, Raymond Smith mailto:smit...@mit.edu>> wrote: Hi, Krishna. Perhaps I'm misunderstanding something, but I'm still not convinced the second version you suggested -- c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue - c_star))], mesh.facesTop) -- isn't working like you want. Could you look at the example I suggested to see if that behaves differently than you expect? Here's the code I used. To me it looks very similar in form to c constraint above and at first glance it seems to behave exactly like we want -- that is, throughout the time stepping, n*grad(phi) is constrained to the value, -phi at the surface. Correct me if I'm wrong, but my impression is that this is the behavior you desire. from fipy import * nx = 50 dx = 1. mesh = Grid1D(nx=nx, dx=dx) phi = CellVariable(name="field variable", mesh=mesh, value=1.0) D = 1. valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is the problema
RE: casting implicit Boundary Conditions in FiPy
Hi Dan, Thank you for the clarifications. Indeed they were quite helpful in getting a more concrete understanding of the issues "I'm confused by your question. Backward Euler refers to the time stepping scheme, which is indeed first order. The spatial discretization scheme is second order and using the boundary condition trick that I demonstrated is second order and also implicit." ** yes. What I am conveying is that, the spatial equation phi_at_boundary = phi_at_lastcell + (gradient/flux at boundary)*(distance between the boundary and last cell-node), is analogous to future_value = present_value + (slope at future value * delta_t), a Backward Euler step in the time-domain. n.grad(phi) = k * phi_at_boundary = k * (phi_at_lastcell + (gradient/flux at boundary)*(distance between the boundary and last cell-node)), might be a first order approximation ? Again, I am unsure of this claim. At a very first glance, I thought the derivation was similar/analgous to Backward Euler Well, now getting into real issue/question. I am trying to extend this concept on a pseudo-2D domain. Let me explain. 1. The diffusion process happens only along the y-axis, ie. we have anisotropic diffusion. 2. The same implicit Neumann BC applies for each x co-ordinate, i.e. we have d_phi/dy along top surface = - k * phi at the top surface My code is shown below. I think that your derivation also holds good in this pseudo-2D diffusion case. With that sorted out, 1. I substituted the 2D equivalent constructions instead of their 1D counterparts 2. For this demo case, I chose a 2x2 grid, with just nx=2 and ny=2 and unit spacing in each direction 3. Defined the diffusion coefficient as a rank 2 faceVariable, and set its value to a suitable tensor, such as to restrict movement in 2D direction. 4. Constrain the diffusion coefficient along the top boundary to be zero. This boundary is where the special implicit Neumann BC applies 5. Calculate the implicitCoeff with the same formula/code (since the derivation is the same) 6. Define the source term coefficient as (mesh.faceNormals * implicitCoeff * mesh.facesTop).divergence) Now, my problem is the following. Although the simulation seems to yield a plausible-looking viewer animation, showing diffusion along y-axis, the mesh.faceNormals * implicitCoeff * mesh.facesTop).divergence return a 'AddoverFaceVariable' object of a much reduced dimension. The mesh.faceNormals * implicitCoeff * mesh.facesTop correctly returns the expected vectors of size 12(for this 2x2 problem, with unit spacing, there are 12 faceCenters). It correctly sets only the top face to the implicit coefficient.But the moment, the divergence operator is called upon, the size of the returned 'AddoverFaceVariable' object drops to 4, with the last 2 values being set to the ImplicitCoeff. Also, it just becomes a 1D vector, with first two elements being zero. Shouldn't this ImplicitSourceTerm coefficient remain as a 2D matrix ? The code runs, but I am worried if I am correctly implementing the boundary condition. Is there any (recommended/standard) way to validate my solution ? I don't have the analytical solution. import fipy as fp nx = 2 ny = nx dx = 1. dy = dx mesh = fp.Grid2D(dx=dx, dy=dy, nx=nx, ny=ny) phi = fp.CellVariable(name = "Field Variable",mesh = mesh) valueBottom = 1. phi.constrain(valueBottom, mesh.facesBottom) D = 1.0 diffCoeff = fp.FaceVariable(mesh=mesh, rank=2, value=[[[0,0],[0,D]]]) diffCoeff.constrain(0., mesh.facesTop) k = -1.0 # We are modelling the BC (d_phi/dx = k*phi) at the top surface implicitCoeff = D*k / (1. - k*dx/2.) # in order to model implicit BC # print (mesh.faceNormals * implicitCoeff * mesh.facesTop).divergence eq = (fp.TransientTerm() == fp.DiffusionTerm(diffCoeff) + fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * mesh.facesTop).divergence)) viewer = fp.Viewer(vars=phi, datamin=0., datamax=1.) viewer.plot() Transient Simulation timeStep = 0.9 * dx**2 / (2 * D) steps = 75 for step in range(steps): eq.solve(var=phi, dt=timeStep) viewer.plot() Krishna -Original Message- From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 13 June 2016 03:17 To: Multiple recipients of list Subject: Re: casting implicit Boundary Conditions in FiPy On Sun, Jun 12, 2016 at 8:14 AM, Gopalakrishnan, Krishnakumar mailto:k.gopalakrishna...@imperial.ac.uk>> wrote: > > 1. The backward Euler method is used in formulating the equation, > "n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) " . But Backward Euler > is a first order method, isn't it ? I am a bit confused about the second > order accuracy statement. I'm confused by your question. Backward Euler refers
Re: casting implicit Boundary Conditions in FiPy
On Sun, Jun 12, 2016 at 8:14 AM, Gopalakrishnan, Krishnakumar wrote: > > 1. The backward Euler method is used in formulating the equation, > "n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) " . But Backward Euler > is a first order method, isn't it ? I am a bit confused about the second > order accuracy statement. I'm confused by your question. Backward Euler refers to the time stepping scheme, which is indeed first order. The spatial discretization scheme is second order and using the boundary condition trick that I demonstrated is second order and also implicit. > 2. In, the equation, "eq = (fp.TransientTerm() == > fp.DiffusionTerm(diffCoeff) - \ > > fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * > mesh.facesRight).divergence)) " , > > from what I understand the 'ImplicitCoeff' is zeroed out in all other > locations other than the right boundary by the boolean method > 'mesh.facesRight' , right ? Exactly, mesh.facesRight is just evaluates to a boolean array. There is also mesh.exteriorFaces which gets you all exterior faces and then you can also do logical operations based on physical location to get more specific masks to capture boundaries or internal regions. > If this is being implemented, do we have the restriction to use fixed dx ? I guess not if you're careful, but it could be a bit fiddly getting the correct distance. > Although dx appears in the Implicit term, it can just be set to the value > of 'dx' of the last node (i.e. length of the last segment in a > variable-mesh sized 1D geometry file).Am I missing something here ? No, but you need to know how to get the correct distance. There are arrays of cell to face distances, but I'd have to think about how to get that working right. > 3. Also, I am a bit confused about the negative signs used in the equation > above.The implicitCoeff is given a negative sign, and the equation's > ImplicitSourceTerm is also given a negative sign. From what I undertand, > this is what we intend to do ? I was just flipping signs to make it work so I'm not sure. Maybe it would be better to have them both positive and it would make more sense that way. > Assign the diffusion coefficient D throughout the domain, using the > faceVariable method. > Manually set it to zero at the right boundary > Then **add** the Diffusion achieved by the div.(D grad( phi)) by using the > ImplicitSourceTerm method, and ensure that this effect gets added only at > the right boundary. You got it. > If this is the workflow, then perhaps we just have to declare the > ImplicitSourceTerm to be positive, and use the addition sign instead of > subtraction sign for the implicitSourceTerm in the equation ? Yup, I think so. > 5. Taking one massive step to get to the answer. Here is a big question > (due to my poor knowledge in the subject) > > I understand that the implicitness enables us to take these large time-steps > due to the affine and A-stable properties of Implicit Methods.But are we > required to do this ? Oh no. You can still take as small a time step as necessary to capture the dynamics at the time scale of interest. It's just that you're not limited by the horrible dx**2 / D time step limit associated with explicit diffusion. Typically it's best to resolve to some dx based time scale (rather than dx**2) that captures an interface being tracked or some feature of interest moving around. >My problem is that, I am solving a coupled system > of PDEs and my other equations, require me to step very very slowly, due to > the inherent stiffness of my system. The PDEs describe processes of varying > time-constants. So, does taking the timestep to be small inhibit us > in any way ? The only down side to small time steps is that their small :-) > Once again, thanks a lot for your solution. Based on the discusssions with > Ray and from this particular method from your posting to the list, I > think I shall be able to solve the problem in 1D. The issues you pointed > out about using this approach in 2D went right over my head. I have to > think about this, when I implement them perhaps in a week or so. Good luck with it. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
Hi Thank you for the solution. This indeed looks like a nifty way to solve the issue at hand. Being a novice in FiPy, and Finite Volume methods in general, I am a bit fuzzy about a couple of things. 1. The backward Euler method is used in formulating the equation, "n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) " . But Backward Euler is a first order method, isn't it ? I am a bit confused about the second order accuracy statement. 2. In, the equation, "eq = (fp.TransientTerm() == fp.DiffusionTerm(diffCoeff) - \ fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * mesh.facesRight).divergence)) " , from what I understand the 'ImplicitCoeff' is zeroed out in all other locations other than the right boundary by the boolean method 'mesh.facesRight' , right ? If this is being implemented, do we have the restriction to use fixed dx ? Although dx appears in the Implicit term, it can just be set to the value of 'dx' of the last node (i.e. length of the last segment in a variable-mesh sized 1D geometry file).Am I missing something here ? 3. Also, I am a bit confused about the negative signs used in the equation above.The implicitCoeff is given a negative sign, and the equation's ImplicitSourceTerm is also given a negative sign. From what I undertand, this is what we intend to do ? * Assign the diffusion coefficient D throughout the domain, using the faceVariable method. * Manually set it to zero at the right boundary * Then **add** the Diffusion achieved by the div.(D grad( phi)) by using the ImplicitSourceTerm method, and ensure that this effect gets added only at the right boundary. If this is the workflow, then perhaps we just have to declare the ImplicitSourceTerm to be positive, and use the addition sign instead of subtraction sign for the implicitSourceTerm in the equation ? 4. Yes, the k needs to be negative for the analytical solution to be stable. My apologies. I should have been more clear in my earlier posts to the list. 5. Taking one massive step to get to the answer. Here is a big question (due to my poor knowledge in the subject) I understand that the implicitness enables us to take these large time-steps due to the affine and A-stable properties of Implicit Methods.But are we required to do this ?My problem is that, I am solving a coupled system of PDEs and my other equations, require me to step very very slowly, due to the inherent stiffness of my system. The PDEs describe processes of varying time-constants. So, does taking the timestep to be small inhibit us in any way ? Once again, thanks a lot for your solution. Based on the discusssions with Ray and from this particular method from your posting to the list, I think I shall be able to solve the problem in 1D. The issues you pointed out about using this approach in 2D went right over my head. I have to think about this, when I implement them perhaps in a week or so. Thanks once again, Krishna From: fipy-boun...@nist.gov on behalf of Daniel Wheeler Sent: Friday, June 10, 2016 6:17 PM To: Multiple recipients of list Subject: Re: casting implicit Boundary Conditions in FiPy Hi Krishina and Ray, Thanks for the interesting discussion. I'm not 100% sure about everything that Krishina is asking for in the latter part of the discussion so I'm just going to address the code that Ray has developed below (my code is below). I think there is a way to handle the right hand side boundary condition both implicitly and with second order accuracy using an ImplicitSourceTerm boundary condition trick. Sort of similar to what is here, http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-fixed-flux-boundary-conditions. Like I said, I think the code below is both second order accurate (for fixed dx) and implicit. Extending this to 2D might raise a few issues. The grid spacing needs to be in the coefficient so, obviously, dx needs to be fixed in both x and y directions. Also, I'm not sure if I'm missing a factor of dx in the source term in 1D as I'm using both a divergence and an ImplicitSourceTerm so there is some question about volumes and face areas in 2D as well. I'm also confused about the signs, I had to flip the sign in front of the source to make it work. It seems to look right though in 1D and you can just take one massive time step to get to the answer. This will only work if k is negative otherwise it's unstable, right? The way I came up with the source was the following n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) and then solve for n.grad(phi) which gives n.grad(phi) = k * phi_P / (1 - dx * k / 2) where phi_P is the value of phi at the cell next to the boundary with the implicit boundary condition. Then we fake the outbound flux to be
Re: casting implicit Boundary Conditions in FiPy
resent an (oversimplified) >> analogous problem when posting to the group, retaining the generality, >> since many other subject experts might have faced similar situations. >> >> >> >> The actual problem I am solving is the solid diffusion PDE (only 1 >> equation) in a Li-ion battery. I am solving this PDE in a pseudo-2D >> domain. i.e. I have defined a Cartesian 2D space, wherein the y-coordinate >> corresponds to the radial direction. The bottom face corresponds to >> particle centres, and the top face corresponding to surface of each >> spherical particle. The x-axis co-ordinate corresponds to particles along >> the width (or thickness) of the positive electrode domain. Diffusion of Li >> is restricted to be within the solid particle (i.e. y-direction only), by >> defining a suitable tensor diffusion coefficient as described in the >> Anisotropic diffusion example and FAQ in FiPy. I have normalised my x and >> y dimensions to have a length of unity. >> >> >> >> Now, the boundary condition along the top face is >> >> >> >> Now, j is non-linear (Butler-Volmer), and I am using a Taylor-expanded >> linear version for this boundary condition. All other field variables >> are assumed as constants. The idea is to set up the infrastructure and >> solve this problem independently, before worrying upon the rubrics of >> setting up the coupled system. In a similar fashion, I have built up and >> solved the solid phase potential PDE (thanks to your help for pointing out >> about the implicit source term). Thus, the idea is to build up the coupled >> P2D Newman model piecemeal. >> >> >> >> The linearised version of my BC’s RHS at a given operating point ()is >> >> >> >> As you can see, the linearised Boundary condition, is cast in terms of >> the field variable, . Hence, we need it in an implicit form >> corresponding to (pseudocode: c.faceGrad.constrain([-( (j_at_c_star - >> partial_j_at_op_point*c_star) + coeff = >> partial_j_at_op_point)],mesh.facesTop) , or something of this >> form/meaning. (just like the very useful ImplicitSourceterm method) >> >> >> >> If I instead apply the c.faceValue method,i.e. using it in setting the >> BC as >> >> >> >> c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue >> - c_star))], mesh.facesTop), then c.faceValue gets immediately >> evaluated at the operating point, c_star, and we are left with 0 >> multiplying the first-order derivative. >> >> >> >> ie. the Boundary conditions becomes, >> >> >> >> Leading to huge loss of accuracy. >> >> >> >> Is there any hope at all in this situation ? J . Cheers and thanks for >> your help thus far. >> >> >> >> >> >> Krishna >> >> >> >> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf >> Of *Raymond Smith >> *Sent:* 09 June 2016 16:06 >> >> *To:* fipy@nist.gov >> *Subject:* Re: casting implicit Boundary Conditions in FiPy >> >> >> >> Hi, Krishna. >> >> Could you give a bit more detail and/or an example about how you know >> it's doing the wrong thing throughout the solution process? In the example >> you sent, the correct solution is the same (c(x, t) = 0) whether you set >> n*grad(phi) to zero or to phi at the boundary, so it's not a good example >> for concluding that it's not behaving as you'd expect. It's helpful here to >> find a situation in which you know that analytical solution to confirm one >> way or the other. For example, you should be able to get the solution to >> the following problem using a Fourier series expansion: >> >> dc/dt = Laplacian(c) >> >> c(t=0) = 1 >> >> x=0: c = 0 >> >> x=1: c - dc/dx = 0 >> >> Ray >> >> >> >> On Thu, Jun 9, 2016 at 10:52 AM, Gopalakrishnan, Krishnakumar < >> k.gopalakrishna...@imperial.ac.uk> wrote: >> >> Hi Ray, >> >> >> >> Thanks for your help. >> >> >> >> But when I apply phi.harmonicFaceValue , it is immediately evaluated to >> a numerical result (a zero vector in this case, since initial value = 0, >> the data type is >> fipy.variables.harmonicCellToFaceVariable._HarmonicCellToFaceVariable', >> i.e. the boundary condition is not remaining implicit. >> >> >> >> The same is the case with the examples.convection.robin example. He
Re: casting implicit Boundary Conditions in FiPy
Okay, I see your dilemma more clearly now; thanks for clarifying and yes, I was missing a few things, so this makes more sense now. Anyway, at this point, we're stepping beyond of my FiPy-fu, so I'll be curious to learn as well. Best, Ray On Fri, Jun 10, 2016 at 7:59 AM, Gopalakrishnan, Krishnakumar < k.gopalakrishna...@imperial.ac.uk> wrote: > Hi Ray, > > > > Thanks a lot for your further inputs. > > > > I realise that perhaps I have not been articulating the problem correctly, > and there is a miscommunication. > > > > In your code, you have the k3 term as constant. More importantly, the > boundary constraint code is set outside the loop , i.e. before beginning > time-stepping. > > > > In this case, a single time-step difference will not make a difference, > i.e . an explicit use of ‘j’ from previous time-step will not be too bad > (we’ll have to think a bit more about stability, due to the time-delay > introduced) > > > > However, my problem is quite different. The linearization of my > (non-linear) RHS of the boundary condition is performed around the > operating point c*. The operating point changes at every time-step. Thus, > the line of code implementing the BC constraint must be inside the loop. > The problem now becomes the following: > > > > 1. The code print c.faceValue return a numerical result. That means > that this expression is evaluated numerically, rather than being left as > part of an implicit exprerssion. The only way a numerical solution can be > obtained is to simply return the existing solution, i.e. c.faceValue > simply uses solution values from the previous time-step. > > 2. After computing the solution, the c* point also needs to be > changed to the new operating point, i.e. c* = c. (You always linearise > around your latest operating point – principle of small-signal perturbation > and linearization) > > 3. Thus, in the next iteration and all subsequent time-steps, > c.faceValue – c* term evaluates to zero again. Thus, we are solving the > PDE only with the DC component of the Taylor series expansion (i.e. > effectively ignoring the first-order derivatives too !). > > > > Note that this sharply contrasts with the problem with the implicit PDE > set-up, i.e. one could have an implicit source term, which is a linear > function of solution variable, and fiPy has a built-in method to handle > this in terms of ImplicitSourceTerm(coeff=…). You simply obtain the > linearised mathematical expression of the source term in paper, and feed it > to the coeff term in the ImplicitSourceTerm method. c.faceValue must not > get evaluated numerically, and must be used implicitly as part of obtaining > the solution at this time-step. > > > > I guess describing the problem in words are falling a bit short, since the > idea is complex enough. I have drawn a simple flow-chart of the workflow > for solving this Implicit BC problem. Can you perhaps kindly take a look > at it ? > https://onedrive.live.com/redir?resid=618ADC32B6C2B152!19186&authkey=!AHrlNkAcvENEZ0M&v=3&ithint=photo%2cjpg > > > > Meanwhile, may I also ask if other Fipy users or developers had to deal > with non-linear Implicit Neumann boundary conditions in their problems? > > > > > > Krishna > > > > *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf > Of *Raymond Smith > *Sent:* 10 June 2016 00:16 > > *To:* fipy@nist.gov > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > Hi, Krisnha. > > I don't know the details about how FiPy handles the boundary conditions, > so I can't confirm whether or not it's solving it with the current > (implicit) or previous (explicit) time step's value when constraining the > gradient. I think that will have to come from one of the devs or someone > more involved than I. However, I did note that in the simple example code > below, I put it in the same form as you have with the extra terms, > n*grad(c) = k1 + k2*(c - k3) > where k1, k2, and k3 are all constants. > It also seems to work fine. It may be treated explicitly, but it does at > least seem to update in time correctly. Is there a significant problem with > it being explicit (using previous time step value) anyway? The overall > order of accuracy should be the same as implicit (current time step), > although I guess stability could be worse? > > In other words, for the first time step, the (c - cstar) term _should_ be > very close to zero (exactly zero when treated explicitly), because c is > initially equal to cstar at the surface. If your time step is large enough > that c at the surface deviates &q
Re: casting implicit Boundary Conditions in FiPy
Hi, Krisnha. I don't know the details about how FiPy handles the boundary conditions, so I can't confirm whether or not it's solving it with the current (implicit) or previous (explicit) time step's value when constraining the gradient. I think that will have to come from one of the devs or someone more involved than I. However, I did note that in the simple example code below, I put it in the same form as you have with the extra terms, n*grad(c) = k1 + k2*(c - k3) where k1, k2, and k3 are all constants. It also seems to work fine. It may be treated explicitly, but it does at least seem to update in time correctly. Is there a significant problem with it being explicit (using previous time step value) anyway? The overall order of accuracy should be the same as implicit (current time step), although I guess stability could be worse? In other words, for the first time step, the (c - cstar) term _should_ be very close to zero (exactly zero when treated explicitly), because c is initially equal to cstar at the surface. If your time step is large enough that c at the surface deviates "significantly" from cstar within a single time step, then your time steps are probably too large anyway. Of course, the (c - cstar) should become non-zero at least by the second time step as long as k1 is finite. from fipy import * nx = 50 dx = 1. mesh = Grid1D(nx=nx, dx=dx) phi = CellVariable(name="field variable", mesh=mesh, value=1.0) D = 1. k1 = -0.1 k2 = 0.1 phistar = 1.0 valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain([-phi.harmonicFaceValue], mesh.facesRight) # This is the problematic BC phi.faceGrad.constrain([k1 + k2*(phistar - phi.harmonicFaceValue)], mesh.facesRight) eq = TransientTerm() == DiffusionTerm(coeff=D) timeStep = 0.9 * dx**2 / (2 * D) steps = 800 viewer = Viewer(vars=phi, datamax=1., datamin=0.) for step in range(steps): eq.solve(var=phi, dt=timeStep) viewer.plot() Cheers, Ray On Thu, Jun 9, 2016 at 6:32 PM, Gopalakrishnan, Krishnakumar < k.gopalakrishna...@imperial.ac.uk> wrote: > Hi Ray, > > > > Thank you for replying. Yes, for this problem solved by your code, it > doesn’t have any other terms inside the [] in > phi.faceGrad.constrain([-phi.harmonicFaceValue], > mesh.facesRight). > > > > What is really confusing me is that, the variable to solve for is phi, > right ? And, we need the solution for this particular time-step within > the loop. But the phi.harmonicFaceValue method performs a numerical > evaluation , i.e. substituting the value of field variable (solution) from > the previous time-step. > > > > If this is indeed the case, there is a problem, because unlike the minimal > example which works correctly, the real problem has the term (c.faceValue > - c_star) . Upon seeing the c.faceValue term, fiPy automatically > substitutes the numerical value of the field variable, ie. boundary value > from the previous time-step. However, this is the same as c_star, which > I methodically verified. Thus, the two terms gets cancelled out, > effectively eliminating the first-order derivative term from the > taylor-series expansion of the non-linear j, thereby leading to > first-order truncation errors. > > > > Is my understanding of how the whole process works correct ? What I am > saying is, there should be a way to cast the (linearised) solution > variable term in the BC, in an implicit format (without substituting the > numerical value from the previous time-step), alike the > ImplicitSourceTerm(coeff=…) style of casting. > > > > > > The solution variable ‘c’ must be obtained by solving the PDE, the > left-boundary no-flux BC , and this right boundary *Implicit* BC, at this > time-step, ie. without substituting c.facevalue (numerical values from > previous time-step) > > > > If we do this c.faceValue substitution, we are simply left with > > > > > > i.e. evaluation of the j at a given DC operating point, introducing first > order errors in the local time-step. > > > > Am I on the wrong track here ? > > > > > > Krishna > > > > *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf > Of *Raymond Smith > *Sent:* 09 June 2016 17:03 > *To:* fipy@nist.gov > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > Hi, Krishna. > > Perhaps I'm misunderstanding something, but I'm still not convinced the > second version you suggested -- c.faceGrad.constrain([-(j_at_c_star + > partial_j_at_op_point*(c.faceValue - c_star))], mesh.facesTop) -- isn't > wo
RE: casting implicit Boundary Conditions in FiPy
Hi Ray, Thank you for replying. Yes, for this problem solved by your code, it doesn’t have any other terms inside the [] in phi.faceGrad.constrain([-phi.harmonicFaceValue], mesh.facesRight). What is really confusing me is that, the variable to solve for is phi,right ? And, we need the solution for this particular time-step within the loop. But the phi.harmonicFaceValue method performs a numerical evaluation , i.e. substituting the value of field variable (solution) from the previous time-step. If this is indeed the case, there is a problem, because unlike the minimal example which works correctly, the real problem has the term (c.faceValue - c_star) . Upon seeing the c.faceValue term, fiPy automatically substitutes the numerical value of the field variable, ie. boundary value from the previous time-step. However, this is the same as c_star, which I methodically verified. Thus, the two terms gets cancelled out, effectively eliminating the first-order derivative term from the taylor-series expansion of the non-linear j, thereby leading to first-order truncation errors. Is my understanding of how the whole process works correct ? What I am saying is, there should be a way to cast the (linearised) solution variable term in the BC, in an implicit format (without substituting the numerical value from the previous time-step), alike the ImplicitSourceTerm(coeff=…) style of casting. [cid:image010.png@01D1C2A7.36988FF0] The solution variable ‘c’ must be obtained by solving the PDE, the left-boundary no-flux BC , and this right boundary Implicit BC, at this time-step, ie. without substituting c.facevalue (numerical values from previous time-step) If we do this c.faceValue substitution, we are simply left with [cid:image011.png@01D1C2A7.36988FF0] i.e. evaluation of the j at a given DC operating point, introducing first order errors in the local time-step. Am I on the wrong track here ? Krishna From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Raymond Smith Sent: 09 June 2016 17:03 To: fipy@nist.gov Subject: Re: casting implicit Boundary Conditions in FiPy Hi, Krishna. Perhaps I'm misunderstanding something, but I'm still not convinced the second version you suggested -- c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue - c_star))], mesh.facesTop) -- isn't working like you want. Could you look at the example I suggested to see if that behaves differently than you expect? Here's the code I used. To me it looks very similar in form to c constraint above and at first glance it seems to behave exactly like we want -- that is, throughout the time stepping, n*grad(phi) is constrained to the value, -phi at the surface. Correct me if I'm wrong, but my impression is that this is the behavior you desire. from fipy import * nx = 50 dx = 1. mesh = Grid1D(nx=nx, dx=dx) phi = CellVariable(name="field variable", mesh=mesh, value=1.0) D = 1. valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is the problematic BC phi.faceGrad.constrain([-phi.harmonicFaceValue], mesh.facesRight) # This is the problematic BC eq = TransientTerm() == DiffusionTerm(coeff=D) timeStep = 0.9 * dx**2 / (2 * D) steps = 800 viewer = Viewer(vars=phi, datamax=1., datamin=0.) for step in range(steps): eq.solve(var=phi, dt=timeStep) viewer.plot() Cheers, Ray On Thu, Jun 9, 2016 at 11:44 AM, Gopalakrishnan, Krishnakumar mailto:k.gopalakrishna...@imperial.ac.uk>> wrote: Hi Ray, Yes. You make a good point. I see that the analytical solution to the particular problem I have posted is also zero. The reason I posted is because I wanted to present an (oversimplified) analogous problem when posting to the group, retaining the generality, since many other subject experts might have faced similar situations. The actual problem I am solving is the solid diffusion PDE (only 1 equation) in a Li-ion battery. I am solving this PDE in a pseudo-2D domain. i.e. I have defined a Cartesian 2D space, wherein the y-coordinate corresponds to the radial direction. The bottom face corresponds to particle centres, and the top face corresponding to surface of each spherical particle. The x-axis co-ordinate corresponds to particles along the width (or thickness) of the positive electrode domain. Diffusion of Li is restricted to be within the solid particle (i.e. y-direction only), by defining a suitable tensor diffusion coefficient as described in the Anisotropic diffusion example and FAQ in FiPy. I have normalised my x and y dimensions to have a length of unity. Now, the boundary condition along the top face is [cid:image001.png@01D1C290.47669780] Now, j is non-linear (Butler-Volmer), and I am using a Taylor-expanded linear version for this bound
Re: casting implicit Boundary Conditions in FiPy
Hi, Krishna. Perhaps I'm misunderstanding something, but I'm still not convinced the second version you suggested -- c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue - c_star))], mesh.facesTop) -- isn't working like you want. Could you look at the example I suggested to see if that behaves differently than you expect? Here's the code I used. To me it looks very similar in form to c constraint above and at first glance it seems to behave exactly like we want -- that is, throughout the time stepping, n*grad(phi) is constrained to the value, -phi at the surface. Correct me if I'm wrong, but my impression is that this is the behavior you desire. from fipy import * nx = 50 dx = 1. mesh = Grid1D(nx=nx, dx=dx) phi = CellVariable(name="field variable", mesh=mesh, value=1.0) D = 1. valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is the problematic BC phi.faceGrad.constrain([-phi.harmonicFaceValue], mesh.facesRight) # This is the problematic BC eq = TransientTerm() == DiffusionTerm(coeff=D) timeStep = 0.9 * dx**2 / (2 * D) steps = 800 viewer = Viewer(vars=phi, datamax=1., datamin=0.) for step in range(steps): eq.solve(var=phi, dt=timeStep) viewer.plot() Cheers, Ray On Thu, Jun 9, 2016 at 11:44 AM, Gopalakrishnan, Krishnakumar < k.gopalakrishna...@imperial.ac.uk> wrote: > Hi Ray, > > > > Yes. You make a good point. I see that the analytical solution to the > particular problem I have posted is also zero. > > > > The reason I posted is because I wanted to present an (oversimplified) > analogous problem when posting to the group, retaining the generality, > since many other subject experts might have faced similar situations. > > > > The actual problem I am solving is the solid diffusion PDE (only 1 > equation) in a Li-ion battery. I am solving this PDE in a pseudo-2D > domain. i.e. I have defined a Cartesian 2D space, wherein the y-coordinate > corresponds to the radial direction. The bottom face corresponds to > particle centres, and the top face corresponding to surface of each > spherical particle. The x-axis co-ordinate corresponds to particles along > the width (or thickness) of the positive electrode domain. Diffusion of Li > is restricted to be within the solid particle (i.e. y-direction only), by > defining a suitable tensor diffusion coefficient as described in the > Anisotropic diffusion example and FAQ in FiPy. I have normalised my x and > y dimensions to have a length of unity. > > > > Now, the boundary condition along the top face is > > > > Now, j is non-linear (Butler-Volmer), and I am using a Taylor-expanded > linear version for this boundary condition. All other field variables are > assumed as constants. The idea is to set up the infrastructure and solve > this problem independently, before worrying upon the rubrics of setting up > the coupled system. In a similar fashion, I have built up and solved the > solid phase potential PDE (thanks to your help for pointing out about the > implicit source term). Thus, the idea is to build up the coupled P2D Newman > model piecemeal. > > > > The linearised version of my BC’s RHS at a given operating point ()is > > > > As you can see, the linearised Boundary condition, is cast in terms of the > field variable, . Hence, we need it in an implicit form corresponding > to (pseudocode: c.faceGrad.constrain([-( (j_at_c_star - > partial_j_at_op_point*c_star) + coeff = > partial_j_at_op_point)],mesh.facesTop) , or something of this > form/meaning. (just like the very useful ImplicitSourceterm method) > > > > If I instead apply the c.faceValue method,i.e. using it in setting the BC > as > > > > c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue - > c_star))], mesh.facesTop), then c.faceValue gets immediately evaluated > at the operating point, c_star, and we are left with 0 multiplying the > first-order derivative. > > > > ie. the Boundary conditions becomes, > > > > Leading to huge loss of accuracy. > > > > Is there any hope at all in this situation ? J . Cheers and thanks for > your help thus far. > > > > > > Krishna > > > > *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf > Of *Raymond Smith > *Sent:* 09 June 2016 16:06 > > *To:* fipy@nist.gov > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > Hi, Krishna. > > Could you give a bit more detail and/or an example about how you know it's > doing the wrong thing throughout the solution process? In
RE: casting implicit Boundary Conditions in FiPy
Hi Ray, Yes. You make a good point. I see that the analytical solution to the particular problem I have posted is also zero. The reason I posted is because I wanted to present an (oversimplified) analogous problem when posting to the group, retaining the generality, since many other subject experts might have faced similar situations. The actual problem I am solving is the solid diffusion PDE (only 1 equation) in a Li-ion battery. I am solving this PDE in a pseudo-2D domain. i.e. I have defined a Cartesian 2D space, wherein the y-coordinate corresponds to the radial direction. The bottom face corresponds to particle centres, and the top face corresponding to surface of each spherical particle. The x-axis co-ordinate corresponds to particles along the width (or thickness) of the positive electrode domain. Diffusion of Li is restricted to be within the solid particle (i.e. y-direction only), by defining a suitable tensor diffusion coefficient as described in the Anisotropic diffusion example and FAQ in FiPy. I have normalised my x and y dimensions to have a length of unity. Now, the boundary condition along the top face is [cid:image005.png@01D1C26E.2386E850] Now, j is non-linear (Butler-Volmer), and I am using a Taylor-expanded linear version for this boundary condition. All other field variables[cid:image006.png@01D1C26E.2386E850] are assumed as constants. The idea is to set up the infrastructure and solve this problem independently, before worrying upon the rubrics of setting up the coupled system. In a similar fashion, I have built up and solved the solid phase potential PDE (thanks to your help for pointing out about the implicit source term). Thus, the idea is to build up the coupled P2D Newman model piecemeal. The linearised version of my BC’s RHS at a given operating point ([cid:image007.png@01D1C26E.2386E850])is [cid:image008.png@01D1C26E.2386E850] As you can see, the linearised Boundary condition, is cast in terms of the field variable, [cid:image007.png@01D1C26E.2386E850] . Hence, we need it in an implicit form corresponding to (pseudocode: c.faceGrad.constrain([-( (j_at_c_star - partial_j_at_op_point*c_star) + coeff = partial_j_at_op_point)],mesh.facesTop) , or something of this form/meaning. (just like the very useful ImplicitSourceterm method) If I instead apply the c.faceValue method,i.e. using it in setting the BC as c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue - c_star))], mesh.facesTop), then c.faceValue gets immediately evaluated at the operating point, c_star, and we are left with 0 multiplying the first-order derivative. ie. the Boundary conditions becomes, [cid:image009.png@01D1C26E.2386E850] Leading to huge loss of accuracy. Is there any hope at all in this situation ? ☺ . Cheers and thanks for your help thus far. Krishna From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Raymond Smith Sent: 09 June 2016 16:06 To: fipy@nist.gov Subject: Re: casting implicit Boundary Conditions in FiPy Hi, Krishna. Could you give a bit more detail and/or an example about how you know it's doing the wrong thing throughout the solution process? In the example you sent, the correct solution is the same (c(x, t) = 0) whether you set n*grad(phi) to zero or to phi at the boundary, so it's not a good example for concluding that it's not behaving as you'd expect. It's helpful here to find a situation in which you know that analytical solution to confirm one way or the other. For example, you should be able to get the solution to the following problem using a Fourier series expansion: dc/dt = Laplacian(c) c(t=0) = 1 x=0: c = 0 x=1: c - dc/dx = 0 Ray On Thu, Jun 9, 2016 at 10:52 AM, Gopalakrishnan, Krishnakumar mailto:k.gopalakrishna...@imperial.ac.uk>> wrote: Hi Ray, Thanks for your help. But when I apply phi.harmonicFaceValue , it is immediately evaluated to a numerical result (a zero vector in this case, since initial value = 0, the data type is fipy.variables.harmonicCellToFaceVariable._HarmonicCellToFaceVariable', i.e. the boundary condition is not remaining implicit. The same is the case with the examples.convection.robin example. Here, the phi.faceValue method is used. However, this also results in an immediate numerical evaluation. However, what is actually required is that, the BC must remain implicit (in variable form, without getting numerically evaluated), being cast in terms of the field variable being solved for. Then the solver needs to solve the PDE on the domain to yield the solution of the field variable. [cid:image010.png@01D1C26E.2386E850] I think we need to solve for the PDE, keeping this implicit BC, rather than immediately evaluating the term [cid:image011.png@01D1C26E.2386E850] , since [cid:image012.png@01D1C26E.2386E850] is the field variable to be solved for, i.e. there ought to be some way to cast the Boun
Re: casting implicit Boundary Conditions in FiPy
Hm, changing the right boundary condition to c + dc/dx = 0 may make a somewhat more reasonable output. Going back to my original question, when I change your script to solve that problem, the solution at least looks _reasonable_ to me, although again, worth checking against an analytical solution. That is, the solution value at the right decreases and the slope at the right seems to decrease proportionally. On Thu, Jun 9, 2016 at 11:06 AM, Raymond Smith wrote: > Hi, Krishna. > > Could you give a bit more detail and/or an example about how you know it's > doing the wrong thing throughout the solution process? In the example you > sent, the correct solution is the same (c(x, t) = 0) whether you set > n*grad(phi) to zero or to phi at the boundary, so it's not a good example > for concluding that it's not behaving as you'd expect. It's helpful here to > find a situation in which you know that analytical solution to confirm one > way or the other. For example, you should be able to get the solution to > the following problem using a Fourier series expansion: > > dc/dt = Laplacian(c) > c(t=0) = 1 > x=0: c = 0 > x=1: c - dc/dx = 0 > > Ray > > > On Thu, Jun 9, 2016 at 10:52 AM, Gopalakrishnan, Krishnakumar < > k.gopalakrishna...@imperial.ac.uk> wrote: > >> Hi Ray, >> >> >> >> Thanks for your help. >> >> >> >> But when I apply phi.harmonicFaceValue , it is immediately evaluated to >> a numerical result (a zero vector in this case, since initial value = 0, >> the data type is >> fipy.variables.harmonicCellToFaceVariable._HarmonicCellToFaceVariable', >> i.e. the boundary condition is not remaining implicit. >> >> >> >> The same is the case with the examples.convection.robin example. Here, >> the phi.faceValue method is used. However, this also results in an >> immediate numerical evaluation. >> >> >> >> However, what is actually required is that, the BC must remain implicit >> (in variable form, without getting numerically evaluated), being cast in >> terms of the field variable being solved for. Then the solver needs to >> solve the PDE on the domain to yield the solution of the field variable. >> >> >> >> >> >> I think we need to solve for the PDE, keeping this implicit BC, rather >> than immediately evaluating the term , since is the field variable to be >> solved for, i.e. there ought to be some way to cast the Boundary condition >> as implicit. >> >> >> >> In FiPy, I have previously set up an implicit source term,by using >> the following code snippet, ImplicitSourceTerm(coeff=k) . Perhaps there >> might be an equivalent method in FiPy to set up the implicit BC, I think ? >> >> >> >> >> >> Krishna >> >> >> >> >> >> >> >> >> >> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf >> Of *Raymond Smith >> *Sent:* 09 June 2016 14:23 >> >> *To:* fipy@nist.gov >> *Subject:* Re: casting implicit Boundary Conditions in FiPy >> >> >> >> Oh, right, the boundary condition is applied on a face, so you need the >> facevalue of phi: >> >> phi.faceGrad.constrain([phi.harmonicFaceValue]) >> >> Ray >> >> >> >> On Thu, Jun 9, 2016 at 7:28 AM, Gopalakrishnan, Krishnakumar < >> k.gopalakrishna...@imperial.ac.uk> wrote: >> >> Hi ray, >> >> >> >> Casting the implicit PDE does not work for my problem. FiPy throws up a >> ton of errors. >> >> I am attaching a minimal example (based off example1.mesh.1D) >> >> >> >> *from *fipy *import ** >> >> nx = 50 >> dx = 1. >> mesh = Grid1D(nx=nx, dx=dx) >> >> phi = CellVariable(name=*"field variable"*, mesh=mesh, value=0.0) >> D = 1. >> >> valueLeft = 0.0 >> phi.constrain(valueLeft, mesh.facesLeft) >> phi.faceGrad.constrain([phi], mesh.facesRight) >> >> *# This is the problematic BC *eq = TransientTerm() == DiffusionTerm( >> coeff=D) >> >> timeStep = 0.9 * dx**2 / (2 * D) >> steps = 100 >> >> viewer = Viewer(vars=phi) >> >> *for *step *in *range(steps): >> eq.solve(var=phi, dt=timeStep) >> viewer.plot() >> >> >> >> The errors are as follows: >> >> >> >> line 22, in >> >> eq.solve(var=phi, dt=timeStep) >> >> \fipy\terms\term.py", line 211, in solve >> >> solver = self._prepare
Re: casting implicit Boundary Conditions in FiPy
Hi, Krishna. Could you give a bit more detail and/or an example about how you know it's doing the wrong thing throughout the solution process? In the example you sent, the correct solution is the same (c(x, t) = 0) whether you set n*grad(phi) to zero or to phi at the boundary, so it's not a good example for concluding that it's not behaving as you'd expect. It's helpful here to find a situation in which you know that analytical solution to confirm one way or the other. For example, you should be able to get the solution to the following problem using a Fourier series expansion: dc/dt = Laplacian(c) c(t=0) = 1 x=0: c = 0 x=1: c - dc/dx = 0 Ray On Thu, Jun 9, 2016 at 10:52 AM, Gopalakrishnan, Krishnakumar < k.gopalakrishna...@imperial.ac.uk> wrote: > Hi Ray, > > > > Thanks for your help. > > > > But when I apply phi.harmonicFaceValue , it is immediately evaluated to a > numerical result (a zero vector in this case, since initial value = 0, the > data type is > fipy.variables.harmonicCellToFaceVariable._HarmonicCellToFaceVariable', > i.e. the boundary condition is not remaining implicit. > > > > The same is the case with the examples.convection.robin example. Here, > the phi.faceValue method is used. However, this also results in an > immediate numerical evaluation. > > > > However, what is actually required is that, the BC must remain implicit > (in variable form, without getting numerically evaluated), being cast in > terms of the field variable being solved for. Then the solver needs to > solve the PDE on the domain to yield the solution of the field variable. > > > > > > I think we need to solve for the PDE, keeping this implicit BC, rather > than immediately evaluating the term , since is the field variable to be > solved for, i.e. there ought to be some way to cast the Boundary condition > as implicit. > > > > In FiPy, I have previously set up an implicit source term,by using > the following code snippet, ImplicitSourceTerm(coeff=k) . Perhaps there > might be an equivalent method in FiPy to set up the implicit BC, I think ? > > > > > > Krishna > > > > > > > > > > *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf > Of *Raymond Smith > *Sent:* 09 June 2016 14:23 > > *To:* fipy@nist.gov > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > Oh, right, the boundary condition is applied on a face, so you need the > facevalue of phi: > > phi.faceGrad.constrain([phi.harmonicFaceValue]) > > Ray > > > > On Thu, Jun 9, 2016 at 7:28 AM, Gopalakrishnan, Krishnakumar < > k.gopalakrishna...@imperial.ac.uk> wrote: > > Hi ray, > > > > Casting the implicit PDE does not work for my problem. FiPy throws up a > ton of errors. > > I am attaching a minimal example (based off example1.mesh.1D) > > > > *from *fipy *import ** > > nx = 50 > dx = 1. > mesh = Grid1D(nx=nx, dx=dx) > > phi = CellVariable(name=*"field variable"*, mesh=mesh, value=0.0) > D = 1. > > valueLeft = 0.0 > phi.constrain(valueLeft, mesh.facesLeft) > phi.faceGrad.constrain([phi], mesh.facesRight) > > *# This is the problematic BC *eq = TransientTerm() == DiffusionTerm(coeff > =D) > > timeStep = 0.9 * dx**2 / (2 * D) > steps = 100 > > viewer = Viewer(vars=phi) > > *for *step *in *range(steps): > eq.solve(var=phi, dt=timeStep) > viewer.plot() > > > > The errors are as follows: > > > > line 22, in > > eq.solve(var=phi, dt=timeStep) > > \fipy\terms\term.py", line 211, in solve > > solver = self._prepareLinearSystem(var, solver, boundaryConditions, > dt) > > \fipy\terms\term.py", line 170, in _prepareLinearSystem > > buildExplicitIfOther=self._buildExplcitIfOther) > > \fipy\terms\binaryTerm.py", line 68, in _buildAndAddMatrices > > buildExplicitIfOther=buildExplicitIfOther) > > \fipy\terms\unaryTerm.py", line 99, in _buildAndAddMatrices > > diffusionGeomCoeff=diffusionGeomCoeff) > > \fipy\terms\abstractDiffusionTerm.py", line 337, in _buildMatrix > > nthCoeffFaceGrad = coeff[numerix.newaxis] * > var.faceGrad[:,numerix.newaxis] > > \fipy\variables\variable.py", line 1575, in __getitem__ > > unit=self.unit, > > \fipy\variables\variable.py", line 255, in _getUnit > > return self._extractUnit(self.value) > > \fipy\variables\variable.py", line 561, in _getValue > > value[..., mask] = numerix.array(constraint.value)[..., mask] > > IndexError: index 50 is out of bounds for axis 1 with size 50 > > > > I
RE: casting implicit Boundary Conditions in FiPy
Hi Ray, Thanks for your help. But when I apply phi.harmonicFaceValue , it is immediately evaluated to a numerical result (a zero vector in this case, since initial value = 0, the data type is fipy.variables.harmonicCellToFaceVariable._HarmonicCellToFaceVariable', i.e. the boundary condition is not remaining implicit. The same is the case with the examples.convection.robin example. Here, the phi.faceValue method is used. However, this also results in an immediate numerical evaluation. However, what is actually required is that, the BC must remain implicit (in variable form, without getting numerically evaluated), being cast in terms of the field variable being solved for. Then the solver needs to solve the PDE on the domain to yield the solution of the field variable. [cid:image001.png@01D1C266.EC46CA10] I think we need to solve for the PDE, keeping this implicit BC, rather than immediately evaluating the term [cid:image002.png@01D1C266.EC46CA10] , since [cid:image003.png@01D1C266.EC46CA10] is the field variable to be solved for, i.e. there ought to be some way to cast the Boundary condition as implicit. In FiPy, I have previously set up an implicit source term, [cid:image004.png@01D1C266.EC46CA10]by using the following code snippet, ImplicitSourceTerm(coeff=k) . Perhaps there might be an equivalent method in FiPy to set up the implicit BC, I think ? Krishna From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Raymond Smith Sent: 09 June 2016 14:23 To: fipy@nist.gov Subject: Re: casting implicit Boundary Conditions in FiPy Oh, right, the boundary condition is applied on a face, so you need the facevalue of phi: phi.faceGrad.constrain([phi.harmonicFaceValue]) Ray On Thu, Jun 9, 2016 at 7:28 AM, Gopalakrishnan, Krishnakumar mailto:k.gopalakrishna...@imperial.ac.uk>> wrote: Hi ray, Casting the implicit PDE does not work for my problem. FiPy throws up a ton of errors. I am attaching a minimal example (based off example1.mesh.1D) from fipy import * nx = 50 dx = 1. mesh = Grid1D(nx=nx, dx=dx) phi = CellVariable(name="field variable", mesh=mesh, value=0.0) D = 1. valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC eq = TransientTerm() == DiffusionTerm(coeff=D) timeStep = 0.9 * dx**2 / (2 * D) steps = 100 viewer = Viewer(vars=phi) for step in range(steps): eq.solve(var=phi, dt=timeStep) viewer.plot() The errors are as follows: line 22, in eq.solve(var=phi, dt=timeStep) \fipy\terms\term.py", line 211, in solve solver = self._prepareLinearSystem(var, solver, boundaryConditions, dt) \fipy\terms\term.py", line 170, in _prepareLinearSystem buildExplicitIfOther=self._buildExplcitIfOther) \fipy\terms\binaryTerm.py", line 68, in _buildAndAddMatrices buildExplicitIfOther=buildExplicitIfOther) \fipy\terms\unaryTerm.py", line 99, in _buildAndAddMatrices diffusionGeomCoeff=diffusionGeomCoeff) \fipy\terms\abstractDiffusionTerm.py", line 337, in _buildMatrix nthCoeffFaceGrad = coeff[numerix.newaxis] * var.faceGrad[:,numerix.newaxis] \fipy\variables\variable.py", line 1575, in __getitem__ unit=self.unit, \fipy\variables\variable.py", line 255, in _getUnit return self._extractUnit(self.value) \fipy\variables\variable.py", line 561, in _getValue value[..., mask] = numerix.array(constraint.value)[..., mask] IndexError: index 50 is out of bounds for axis 1 with size 50 I have tried including the implicit BC within the time-stepper loop, but that does not still help. Best Regards Krishna From: fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov> [mailto:fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov>] On Behalf Of Gopalakrishnan, Krishnakumar Sent: 08 June 2016 23:42 To: fipy@nist.gov<mailto:fipy@nist.gov> Subject: RE: casting implicit Boundary Conditions in FiPy Hi Raymond, Sorry, it was a typo. Yes, It is indeed d (phi)/dx, the spatial derivative BC. I shall try setting phi.faceGrad.constrain([k*phi], mesh.facesRight), and see if it will work. Thanks for pointing this out. Krishna From: fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov> [mailto:fipy-boun...@nist.gov] On Behalf Of Raymond Smith Sent: 08 June 2016 23:36 To: fipy@nist.gov<mailto:fipy@nist.gov> Subject: Re: casting implicit Boundary Conditions in FiPy Hi, Krishna. Just to make sure, do you mean that the boundary condition is a derivative with respect to the spatial variable or with respect to time as-written? If you mean spatial, such that d\phi/dx = k*phi, have you tried phi.faceGrad.constrain(k*phi) and that didn't work? If you mean that its value is prescribed by its rate of change, then I'm not sure the best way to do it. Could you maybe do it explicitly? - Store the values from the last time step with hasOld set to True in the crea
Re: casting implicit Boundary Conditions in FiPy
Oh, right, the boundary condition is applied on a face, so you need the
facevalue of phi:
phi.faceGrad.constrain([phi.harmonicFaceValue])
Ray
On Thu, Jun 9, 2016 at 7:28 AM, Gopalakrishnan, Krishnakumar <
k.gopalakrishna...@imperial.ac.uk> wrote:
> Hi ray,
>
>
>
> Casting the implicit PDE does not work for my problem. FiPy throws up a
> ton of errors.
>
> I am attaching a minimal example (based off example1.mesh.1D)
>
>
>
> *from *fipy *import **
>
> nx = 50
> dx = 1.
> mesh = Grid1D(nx=nx, dx=dx)
>
> phi = CellVariable(name=*"field variable"*, mesh=mesh, value=0.0)
> D = 1.
>
> valueLeft = 0.0
> phi.constrain(valueLeft, mesh.facesLeft)
> phi.faceGrad.constrain([phi], mesh.facesRight)
>
> *# This is the problematic BC *eq = TransientTerm() == DiffusionTerm(coeff
> =D)
>
> timeStep = 0.9 * dx**2 / (2 * D)
> steps = 100
>
> viewer = Viewer(vars=phi)
>
> *for *step *in *range(steps):
> eq.solve(var=phi, dt=timeStep)
> viewer.plot()
>
>
>
> The errors are as follows:
>
>
>
> line 22, in
>
> eq.solve(var=phi, dt=timeStep)
>
> \fipy\terms\term.py", line 211, in solve
>
> solver = self._prepareLinearSystem(var, solver, boundaryConditions,
> dt)
>
> \fipy\terms\term.py", line 170, in _prepareLinearSystem
>
> buildExplicitIfOther=self._buildExplcitIfOther)
>
> \fipy\terms\binaryTerm.py", line 68, in _buildAndAddMatrices
>
> buildExplicitIfOther=buildExplicitIfOther)
>
> \fipy\terms\unaryTerm.py", line 99, in _buildAndAddMatrices
>
> diffusionGeomCoeff=diffusionGeomCoeff)
>
> \fipy\terms\abstractDiffusionTerm.py", line 337, in _buildMatrix
>
> nthCoeffFaceGrad = coeff[numerix.newaxis] *
> var.faceGrad[:,numerix.newaxis]
>
> \fipy\variables\variable.py", line 1575, in __getitem__
>
> unit=self.unit,
>
> \fipy\variables\variable.py", line 255, in _getUnit
>
> return self._extractUnit(self.value)
>
> \fipy\variables\variable.py", line 561, in _getValue
>
> value[..., mask] = numerix.array(constraint.value)[..., mask]
>
> IndexError: index 50 is out of bounds for axis 1 with size 50
>
>
>
> I have tried including the implicit BC within the time-stepper loop, but
> that does not still help.
>
>
>
>
>
> Best Regards
>
>
>
> Krishna
>
>
>
>
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Gopalakrishnan, Krishnakumar
> *Sent:* 08 June 2016 23:42
> *To:* fipy@nist.gov
> *Subject:* RE: casting implicit Boundary Conditions in FiPy
>
>
>
> Hi Raymond,
>
>
>
> Sorry, it was a typo.
>
>
>
> Yes, It is indeed d (phi)/dx, the spatial derivative BC. I shall try
> setting phi.faceGrad.constrain([k*phi], mesh.facesRight), and see if it
> will work.
>
>
>
> Thanks for pointing this out.
>
>
>
>
>
> Krishna
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov
> ] *On Behalf Of *Raymond Smith
> *Sent:* 08 June 2016 23:36
> *To:* fipy@nist.gov
> *Subject:* Re: casting implicit Boundary Conditions in FiPy
>
>
>
> Hi, Krishna.
>
> Just to make sure, do you mean that the boundary condition is a derivative
> with respect to the spatial variable or with respect to time as-written? If
> you mean spatial, such that d\phi/dx = k*phi, have you tried
>
> phi.faceGrad.constrain(k*phi) and that didn't work?
>
> If you mean that its value is prescribed by its rate of change, then I'm
> not sure the best way to do it. Could you maybe do it explicitly?
> - Store the values from the last time step with hasOld set to True in the
> creation of the cell variable
> - In each time step, calculate the backward-Euler time derivative
> manually and then set the value of phi with the phi.constrain method.
>
>
>
> Ray
>
>
>
> On Wed, Jun 8, 2016 at 6:26 PM, Gopalakrishnan, Krishnakumar <
> k.gopalakrishna...@imperial.ac.uk> wrote:
>
> I am trying to solve the standard fickean diffusion equation on a 1D
> uniform mesh in (0,1)
>
>
>
> $$\frac{\partial \phi}{\partial t} = \nabla.(D \nabla\phi)$$
>
>
>
> with a suitable initial value for $\phi(x,t)$.
>
>
>
> The problem is that, one of my boundary conditions is *implicit*, i.e. is
> a function of the field variable being solved for.
>
>
>
> $ \frac{\partial\phi}{\partial t} = k \phi $ , at the right boundary edge,
> k = constant
>
>
>
> The left BC is not a problem, it is just a standard no-flux BC.
>
>
>
> How do I cast this *implicit BC* in FiPy ? Any help/poin
RE: casting implicit Boundary Conditions in FiPy
Hi ray,
Casting the implicit PDE does not work for my problem. FiPy throws up a ton of
errors.
I am attaching a minimal example (based off example1.mesh.1D)
from fipy import *
nx = 50
dx = 1.
mesh = Grid1D(nx=nx, dx=dx)
phi = CellVariable(name="field variable", mesh=mesh, value=0.0)
D = 1.
valueLeft = 0.0
phi.constrain(valueLeft, mesh.facesLeft)
phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC
eq = TransientTerm() == DiffusionTerm(coeff=D)
timeStep = 0.9 * dx**2 / (2 * D)
steps = 100
viewer = Viewer(vars=phi)
for step in range(steps):
eq.solve(var=phi, dt=timeStep)
viewer.plot()
The errors are as follows:
line 22, in
eq.solve(var=phi, dt=timeStep)
\fipy\terms\term.py", line 211, in solve
solver = self._prepareLinearSystem(var, solver, boundaryConditions, dt)
\fipy\terms\term.py", line 170, in _prepareLinearSystem
buildExplicitIfOther=self._buildExplcitIfOther)
\fipy\terms\binaryTerm.py", line 68, in _buildAndAddMatrices
buildExplicitIfOther=buildExplicitIfOther)
\fipy\terms\unaryTerm.py", line 99, in _buildAndAddMatrices
diffusionGeomCoeff=diffusionGeomCoeff)
\fipy\terms\abstractDiffusionTerm.py", line 337, in _buildMatrix
nthCoeffFaceGrad = coeff[numerix.newaxis] * var.faceGrad[:,numerix.newaxis]
\fipy\variables\variable.py", line 1575, in __getitem__
unit=self.unit,
\fipy\variables\variable.py", line 255, in _getUnit
return self._extractUnit(self.value)
\fipy\variables\variable.py", line 561, in _getValue
value[..., mask] = numerix.array(constraint.value)[..., mask]
IndexError: index 50 is out of bounds for axis 1 with size 50
I have tried including the implicit BC within the time-stepper loop, but that
does not still help.
Best Regards
Krishna
From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of
Gopalakrishnan, Krishnakumar
Sent: 08 June 2016 23:42
To: fipy@nist.gov
Subject: RE: casting implicit Boundary Conditions in FiPy
Hi Raymond,
Sorry, it was a typo.
Yes, It is indeed d (phi)/dx, the spatial derivative BC. I shall try setting
phi.faceGrad.constrain([k*phi], mesh.facesRight), and see if it will work.
Thanks for pointing this out.
Krishna
From: fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov>
[mailto:fipy-boun...@nist.gov] On Behalf Of Raymond Smith
Sent: 08 June 2016 23:36
To: fipy@nist.gov<mailto:fipy@nist.gov>
Subject: Re: casting implicit Boundary Conditions in FiPy
Hi, Krishna.
Just to make sure, do you mean that the boundary condition is a derivative with
respect to the spatial variable or with respect to time as-written? If you mean
spatial, such that d\phi/dx = k*phi, have you tried
phi.faceGrad.constrain(k*phi) and that didn't work?
If you mean that its value is prescribed by its rate of change, then I'm not
sure the best way to do it. Could you maybe do it explicitly?
- Store the values from the last time step with hasOld set to True in the
creation of the cell variable
- In each time step, calculate the backward-Euler time derivative manually and
then set the value of phi with the phi.constrain method.
Ray
On Wed, Jun 8, 2016 at 6:26 PM, Gopalakrishnan, Krishnakumar
mailto:k.gopalakrishna...@imperial.ac.uk>>
wrote:
I am trying to solve the standard fickean diffusion equation on a 1D uniform
mesh in (0,1)
$$\frac{\partial \phi}{\partial t} = \nabla.(D \nabla\phi)$$
with a suitable initial value for $\phi(x,t)$.
The problem is that, one of my boundary conditions is implicit, i.e. is a
function of the field variable being solved for.
$ \frac{\partial\phi}{\partial t} = k \phi $ , at the right boundary edge, k =
constant
The left BC is not a problem, it is just a standard no-flux BC.
How do I cast this implicit BC in FiPy ? Any help/pointers will be much
appreciated.
Best regards
Krishna
Imperial College London
___
fipy mailing list
fipy@nist.gov<mailto: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: casting implicit Boundary Conditions in FiPy
Hi Raymond,
Sorry, it was a typo.
Yes, It is indeed d (phi)/dx, the spatial derivative BC. I shall try setting
phi.faceGrad.constrain([k*phi], mesh.facesRight), and see if it will work.
Thanks for pointing this out.
Krishna
From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Raymond
Smith
Sent: 08 June 2016 23:36
To: fipy@nist.gov
Subject: Re: casting implicit Boundary Conditions in FiPy
Hi, Krishna.
Just to make sure, do you mean that the boundary condition is a derivative with
respect to the spatial variable or with respect to time as-written? If you mean
spatial, such that d\phi/dx = k*phi, have you tried
phi.faceGrad.constrain(k*phi) and that didn't work?
If you mean that its value is prescribed by its rate of change, then I'm not
sure the best way to do it. Could you maybe do it explicitly?
- Store the values from the last time step with hasOld set to True in the
creation of the cell variable
- In each time step, calculate the backward-Euler time derivative manually and
then set the value of phi with the phi.constrain method.
Ray
On Wed, Jun 8, 2016 at 6:26 PM, Gopalakrishnan, Krishnakumar
mailto:k.gopalakrishna...@imperial.ac.uk>>
wrote:
I am trying to solve the standard fickean diffusion equation on a 1D uniform
mesh in (0,1)
$$\frac{\partial \phi}{\partial t} = \nabla.(D \nabla\phi)$$
with a suitable initial value for $\phi(x,t)$.
The problem is that, one of my boundary conditions is implicit, i.e. is a
function of the field variable being solved for.
$ \frac{\partial\phi}{\partial t} = k \phi $ , at the right boundary edge, k =
constant
The left BC is not a problem, it is just a standard no-flux BC.
How do I cast this implicit BC in FiPy ? Any help/pointers will be much
appreciated.
Best regards
Krishna
Imperial College London
___
fipy mailing list
fipy@nist.gov<mailto: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: casting implicit Boundary Conditions in FiPy
Hi, Krishna.
Just to make sure, do you mean that the boundary condition is a derivative
with respect to the spatial variable or with respect to time as-written? If
you mean spatial, such that d\phi/dx = k*phi, have you tried
phi.faceGrad.constrain(k*phi) and that didn't work?
If you mean that its value is prescribed by its rate of change, then I'm
not sure the best way to do it. Could you maybe do it explicitly?
- Store the values from the last time step with hasOld set to True in the
creation of the cell variable
- In each time step, calculate the backward-Euler time derivative manually
and then set the value of phi with the phi.constrain method.
Ray
On Wed, Jun 8, 2016 at 6:26 PM, Gopalakrishnan, Krishnakumar <
k.gopalakrishna...@imperial.ac.uk> wrote:
> I am trying to solve the standard fickean diffusion equation on a 1D
> uniform mesh in (0,1)
>
>
>
> $$\frac{\partial \phi}{\partial t} = \nabla.(D \nabla\phi)$$
>
>
>
> with a suitable initial value for $\phi(x,t)$.
>
>
>
> The problem is that, one of my boundary conditions is *implicit*, i.e. is
> a function of the field variable being solved for.
>
>
>
> $ \frac{\partial\phi}{\partial t} = k \phi $ , at the right boundary edge,
> k = constant
>
>
>
> The left BC is not a problem, it is just a standard no-flux BC.
>
>
>
> How do I cast this *implicit BC* in FiPy ? Any help/pointers will be much
> appreciated.
>
>
>
>
>
> Best regards
>
>
>
> Krishna
>
> Imperial College London
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
casting implicit Boundary Conditions in FiPy
I am trying to solve the standard fickean diffusion equation on a 1D uniform
mesh in (0,1)
$$\frac{\partial \phi}{\partial t} = \nabla.(D \nabla\phi)$$
with a suitable initial value for $\phi(x,t)$.
The problem is that, one of my boundary conditions is implicit, i.e. is a
function of the field variable being solved for.
$ \frac{\partial\phi}{\partial t} = k \phi $ , at the right boundary edge, k =
constant
The left BC is not a problem, it is just a standard no-flux BC.
How do I cast this implicit BC in FiPy ? Any help/pointers will be much
appreciated.
Best regards
Krishna
Imperial College London
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
