Re: Diffusion-Convection-Reactive Chemisty

2016-07-21 Thread Daniel DeSantis
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

2016-07-21 Thread Campbell, Ian
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

2016-07-21 Thread Daniel Wheeler
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

2016-07-21 Thread Daniel Wheeler
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 ]