Re: Diffusion-Convection-Reactive Chemisty
So, these changes all helped a lot, both in getting results and my general understanding. Sorry to come back with more questions but I have been trying to extend the model a bit. I've added a second, competing reaction and have tried to incorporate a temperature profile. So far, I can get answers that seem to match what I expect but I have to use strange values to do so... The two problems I have are: 1) The temperature profile requires a k_eff of 5e6 or 5e7 to not increase ridiculously high. 2) The diffusion coefficient has to be relatively high (D=100) in order to get a decent concentration profile. In most models with a reaction that has a decent velocity, you can ignore diffusion all together. That does not seem to be the case here. Perhaps my understanding of how to code these systems is incorrect? In this particular case, is there something I should be doing to make the diffusion term second order? I've included codes and an equation sheet to try and clarify things. As it stands, the coefficients (D, k1, k2, keff) in the code give me a concentration profile I expect but not a temperature one. Thank you, Daniel DeSantis On Thu, Jul 14, 2016 at 10:36 AM, Daniel DeSantis wrote: > Thank you everyone. This helps a lot! > > On Wed, Jul 13, 2016 at 5:04 PM, Guyer, Jonathan E. Dr. (Fed) < > jonathan.gu...@nist.gov> wrote: > >> That should be OK. FiPy automatically maps the constraint onto the >> faceValue of a CellVariable. >> >> > On Jul 13, 2016, at 3:26 PM, Keller, Trevor (Fed) < >> trevor.kel...@nist.gov> wrote: >> > >> > Is the definition of C_a_BC correct? For a 1D grid, is the behavior of >> > C_a.constrain(C_a_BC, where=mesh.facesRight) >> > with a CellVariable instead of a scalar value >> > C_a_BC= C_a_0*(1-X) >> > meaningful? >> > >> > Trevor >> > From: fipy-boun...@nist.gov on behalf of >> Daniel DeSantis >> > Sent: Wednesday, July 13, 2016 4:08:50 PM >> > To: FIPY >> > Subject: Re: Diffusion-Convection-Reactive Chemisty >> > >> > I'm sorry. I was trying to fix the problem, and forgot to put a line >> back in, which was masked by me not clearing a variable value for V. My >> apologies. Please try this code instead. It has the traceback I mentioned >> before. >> > >> > >> > >> > Traceback (most recent call last): >> > >> > File "", line 1, in >> > runfile('C:/Users/ddesantis/Dropbox/PythonScripts/CFD >> Models/ConversionModel.py', >> wdir='C:/Users/ddesantis/Dropbox/PythonScripts/CFD Models') >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", >> line 699, in runfile >> > execfile(filename, namespace) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", >> line 74, in execfile >> > exec(compile(scripttext, filename, 'exec'), glob, loc) >> > >> > File "C:/Users/ddesantis/Dropbox/PythonScripts/CFD >> Models/ConversionModel.py", line 110, in >> > res = Eq.sweep(dt=dt, solver=solver) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\term.py", >> line 236, in sweep >> > solver = self._prepareLinearSystem(var=var, solver=solver, >> boundaryConditions=boundaryConditions, dt=dt) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\term.py", >> line 170, in _prepareLinearSystem >> > buildExplicitIfOther=self._buildExplcitIfOther) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\coupledBinaryTerm.py", >> line 122, in _buildAndAddMatrices >> > buildExplicitIfOther=buildExplicitIfOther) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py", >> line 68, in _buildAndAddMatrices >> > buildExplicitIfOther=buildExplicitIfOther) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py", >> line 68, in _buildAndAddMatrices >> > buildExplicitIfOther=buildExplicitIfOther) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py", >> line 68, in _buildAndAddMatrices >> > buildExplicitIfOther=buildExplicitIfOther) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\unaryTerm.py", >> line 99, in _buildAndAddMatrices >> > diffusionGeomCoeff=diffusionGeomCoeff) >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\abstractConvectionTerm.py", >> line 211, in _buildMatrix >> > self.constraintB = -((1 - alpha) * var.arithmeticFaceValue * >> constraintMask * exteriorCoeff).divergence * mesh.cellVolumes >> > >> > File >> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py", >> line 1151, in __mul__
Complex PDE Coupling; Extending to a Further PDE
Dear Jonathan, Thank you very much for your gist* explanation of how to couple the PDEs along the upper boundary of a 2D mesh. Your source code is invaluable and has helped us to get back on track. We have two follow-up questions as a result: Question 1: As an extension of this problem, we're interested in adding a third PDE (PDE3) to the system (this PDE is similar in form to a standard advection-diffusion equation). This third PDE is defined only along the x-axis at the top of the Cartesian grid, just as with PDE1. A component in the source term of PDE3 is the cellvariable in PDE1. Furthermore, another component in the source term of PDE3 exists in the boundary condition of PDE2. That is to say, PDE1 and PDE3 interact as fully coupled equations along the x-axis at the top of the mesh. PDE1 and PDE3 in-turn interact with PDE2 (anisotropic diffusion) along the top of the mesh only. Our question is as follows; given that you have proposed the uncoupled method with repeated sweeping for solving the system of PDE1 and PDE2, does this preclude us from using the fully-coupled "&" method to couple between PDE1 and PDE3? Phrased differently, if we use the uncoupled sweeping method for PDE1 and PDE2, must we continue to use this approach when introducing more axial-domain PDEs? Question 2: In the FiPy coupled-diffusion example (hyperlink below**), the following statement is made: "The uncoupled method still works, but it can be advantageous to solve the two equations simultaneously." Can you please elaborate on this so that we may better understand what the disadvantages of using the uncoupled, repeated sweep method may be, relative to coupling using the "&" operator? Especially given that in our problem, the full-2D implementation has interactions only at the domain's top-edge and it seems to be wasteful of computing resources? Sincerely, - Ian & Krishna * https://gist.github.com/guyer/bb199559c00f6047d466daa18554d83d ** http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.coupled.html Ian Campbell | PhD Candidate Electrochemical Science & Engineering Group Room 506, City & Guilds Building, Imperial College London, SW7 2AZ, United Kingdom Phone: +44 (0)7449 815 520 | E-mail: i.campbel...@imperial.ac.uk ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Improving non-orthogonality issue in FiPy
Hi Michael, I think this paper explains the issue better than I can, https://www.researchgate.net/publication/263972969_ANALYSIS_OF_THE_NON-ORTHOGONALITY_CORRECTION_OF_FINITE_VOLUME_DISCRETIZATION_ON_UNSTRUCTURED_MESHES Figure 2 shows the control volume and Figure 3 shows the skewness issue. Cheers, Daniel On Wed, Jul 20, 2016 at 4:30 PM, Michael Waters wrote: > Hi, I'm not sure I understand the issue. Maybe you could explain more? I > made a diagram to help. > > Cheers, > > -Mike Waters > > > > On 07/20/2016 02:51 PM, Daniel Wheeler wrote: >> >> FiPy has issues with non-orthogonal meshes (i.e. when the vector >> between the cell center isn't parallel to the face normal). We did >> make an attempt at fixing this, which resulted in two diffusion terms, >> >> >> https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermCorrection.py >> >> and >> >> >> https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermNoCorrection.py >> >> As I recall, the "corrected" diffusion term was quite unstable for >> moderate to high non-orthogonality. >> >> I'd like to address this issue again and I was wondering if anyone had >> opinions regarding the latest / best / most tractable approach for >> dealing with non-orthogonality in CC-FV. I haven't followed the >> literature on this for quite some time. >> > > > ___ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Accuracy of DiffusionTerm for non-uniform mesh
On Wed, Jul 20, 2016 at 5:45 PM, Raymond Smith wrote: > That makes sense. So, here > https://gist.github.com/raybsmith/b0b6ee7c90efdcc35d6a0658319f1a01 > I've changed it so that the error is calculated as > sum( (\phi-phi*)^2 * mesh.cellVolumes ) > and depending on the value I choose for the ratio of dx_{i+1} / dx_i, I get > convergence for exponential spacing ranging from 2nd order (at that ratio = > 1) and decreasing order as that ratio decreases from unity. Thanks for that. I'm not entirely sure, but the integral as written on line 44, https://gist.github.com/raybsmith/b0b6ee7c90efdcc35d6a0658319f1a01#file-fipy_accuracy-py-L44 may still only be first order accurate for non-uniform grids. I'll try and investigate this though and see if it matters. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]