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' <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' <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 <fipy@nist.gov<mailto:fipy@nist.gov>> Subject: Re: casting implicit Boundary Conditions in FiPy On Fri, Jun 17, 2016 at 5:36 PM, Gopalakrishnan, Krishnakumar <k.gopalakrishna...@imperial.ac.uk<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'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<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 ]
