Interesting (perhaps unique) approach to parallelising FiPy sweeps

[Gopalakrishnan, Krishnakumar]

Hello,

In our present system of equations, we sweep 8 loosely coupled PDEs defined in 
5 meshes, until their residues die down to an acceptable value.

Upon profiling the code, it's clear that the sweep function is a major 
computational bottleneck. We tried MPI parallelism, but that didn't work out. 
But to side-step this, recognising the fact that the sweep() functions are 
computed one-by-one in serial, we propose an interesting work-around for an 
alternative style of parallelism.

Instead of dividing the control volumes among multiple CPUs, we instead propose 
to split the 8 different sweep calls (in the sweep loop within a time-step) to 
worker cores. Each sweep can then occur independently (instead of the python 
interpreter waiting for one sweep to finish before beginning another).  I 
suspect that this would be beneficial for those problems with 6-8 sweep calls, 
wherein the sweeps are CPU intensive. Thus, the cost of communication overhead 
between the master process to the child workers is offset by the benefits 
gained by independent sweeping.

>From my fairly underwhelming numerical computing knowledge, this looks like 
>analogous to a Jacobi-style implementation, rather than a Gauss-Seidel style 
>iteration, wherein the updated cell-variable from the sweep is immediately 
>available for all subsequent PDEs of the coupled system within the sweep-loop. 
> In the case of multi-processor sweep, the values are updated using values 
>from last sweep iteration.  Thus, this might be slower (but more stable) to 
>converge, isn't it ?  There is also a penalty incurred by the communication 
>overhead of opening and killing processes.

The questions are : 1) Has anybody tried this before ? 2) Does this sound 
problematic from a technical or practical perspective. 3) Does this sound 
remotely useful.

I have successfully implemented a multiprocessor parallelisation for 
examples.diffusion.coupled (shown below) . The values of the simulation results 
(for v0 & v1) after 100 time-steps are very very close to that of the serial 
sweep. Although this code is several times slower than the serial code, we 
tried this more as a concept demonstrator before embarking to spend significant 
effort to convert our serial code to multiprocessor approach.  If you see a red 
flag, it would be much appreciated if you can help us by pointing it out.

from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer
from multiprocessing import Process, Queue, cpu_count

m = Grid1D(nx=100, Lx=1.)

v0 = CellVariable(mesh=m, hasOld=True, value=0.5)
v1 = CellVariable(mesh=m, hasOld=True, value=0.5)

def boundaryConditions(m,v0,v1):
v0.constrain(0, m.facesLeft)
v0.constrain(1, m.facesRight)
v1.constrain(1, m.facesLeft)
v1.constrain(0, m.facesRight)

## Multiprocessor implementation 
#
def get_res0(res0_queue,v0,v1):
boundaryConditions(m, v0, v1)
eq0 = TransientTerm() == DiffusionTerm(coeff=0.01) - v1.faceGrad.divergence
res0 = eq0.sweep(var=v0, dt=1e-5)
res0_queue.put([res0,v0])
def get_res1(res1_queue,v0,v1):
boundaryConditions(m, v0, v1)
eq1 = TransientTerm() == v0.faceGrad.divergence + DiffusionTerm(coeff=0.01)
res1 = eq1.sweep(var=v1, dt=1e-5)
res1_queue.put([res1,v1])

if __name__ == '__main__':
print "You have " + str(cpu_count()) + " CPU cores in your machine"
print "Simulation end time is 99 seconds.\n"
res0_queue = Queue()
res1_queue = Queue()
for t in range(100):
print "Computing the solution at time, t = " + str(t) + " sec"
v0.updateOld()
v1.updateOld()
res0 = res1 = 1e6
while max(res0,res1) > 0.01:
process_a = Process(target=get_res0, args=(res0_queue,v0,v1))
process_b = Process(target=get_res1, args=(res1_queue,v0,v1))
process_a.start()
process_b.start()
process_a.join()
process_b.join()
while not res0_queue.empty():
res0, v0 =  res0_queue.get()
while not res1_queue.empty():
res1, v1 =  res1_queue.get()
print v0, v1
print res0, res1





___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

RE: Applying Two Different Diffusion Coefficients at Single FV Face

[Gopalakrishnan, Krishnakumar]

Thank you, Dr Guyer and Raymond,

That was quite helpful


Krishna & Ian 

-Original Message-
From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Guyer, 
Jonathan E. Dr. (Fed)
Sent: Friday, September 16, 2016 8:11 PM
To: FIPY 
Subject: Re: Applying Two Different Diffusion Coefficients at Single FV Face

Yes!

> On Sep 16, 2016, at 2:18 PM, Raymond Smith  wrote:
> 
> A side note, here it may actually be more convenient to think about the D in 
> terms of the volumes, so you could also define it as a cell variable with 
> values specified as D1 when x<=1 and D2 when x>1, then use the 
> harmonicFaceValue attribute when you put it into the governing equation to 
> have FiPy do this calculation for you and avoid awkward use of the mesh 
> spacing in the specification of the variable values.
> 
> On Fri, Sep 16, 2016 at 11:12 AM, Raymond Smith  wrote:
> Hi, Ian.
> 
> I don't think there is such a thing as having two different flux 
> coefficients at the same face for the same governing PDE. The flux 
> through a given face is calculated by the coefficient at that face times some 
> approximation of the gradient in a field variable at that face, like D * 
> grad(c), both evaluated at the face, so there is only one coefficient in that 
> expression. And in the FiPy FV approach D doesn't come into play at all 
> except at the faces, so it can't change immediately to the right or left of a 
> face but only from one face to the next.
> 
> That said, perhaps you could try following a common approach used when a FV 
> interface is directly between two bulk regions with different continuum 
> properties -- use some sort of a mean of the continuum transport coefficient 
> in the adjacent volumes as the approximation for the coefficient at the face. 
> I would suggest the harmonic mean (one simple reason: it preserves zero flux 
> through the face if either of the adjacent cells has a zero coefficient). You 
> could try:
> Dint = 2*D1*D2/(D1+D2)
> Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) 
> & (1D_mesh.faceCenters[0] < 1.0+dx/2))) If you do that, you should see a 
> change in the slope of the field variable right near the face between the two 
> bulk regions, which seems to be what you're missing. Of course, ideally your 
> mesh is fine enough that the dx is insignificant in terms of the solution, 
> but this is a common approach to get a slightly more reasonable estimate for 
> what happens between two regions.
> 
> Best,
> Ray
> 
> On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian  
> wrote:
> Hello All,
> 
>  
> 
> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard elliptic 
> diffusion equation (\nabla. (D \nabla \phi) = f is being solved.
> 
>  
> 
> The diffusion coefficient, D, changes sharply at x = 1.0. Although this is 
> quite similar to examples.diffusion.mesh1D,  there is a key difference:
> 
>  
> 
> We are seeking an instantaneous change in the coefficient value at that face, 
> i.e. something of the effect that on one side of the face (i.e. from x = 0 up 
> to and including x = 1.0) diffusion coefficient D has a value D1.  Then there 
> is an abrupt change in diffusion coefficient from x >= 1.0 (very notable, 
> including the face at x = 1.0) and for the remainder of the mesh, diffusion 
> coefficient D has a value D2.  Currently, we have to either decide on whether 
> the single face at x = 1.0 gets the value D1 or D2. However, when we run with 
> either of these cases, the results are not accurate (when compared against a 
> benchmark from a commercial FV solver). We are getting a smooth interpolation 
> across the face , rather than the sharp change correctly predicted by the 
> benchmarking code.
> 
>  
> 
> With the code below, we have only been able to set the coefficient at that 
> middle face to one of either D1 or D2. This doesn’t result in the desired 
> instantaneous change in the coefficient value, but instead, the new 
> coefficient value only applies from the NEXT face (which is a distance dx 
> away) from this critical interior boundary.
> 
>  
> 
> Coefficient = FaceVariable(mesh = 1D_mesh)
>   
> # Define facevariable on mesh
> 
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0))   
>  
> # Set coefficient values in 1st half of mesh
> 
> Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) & 
> (1D_mesh.faceCenters[0] <= 2.0)))# Set coefficient values 
> in 2nd half of mesh
> 
>  
> 
> An alternative with the inequalities adjusted, but that still only permits a 
> single coefficient at that face:
> 
> Coefficient = FaceVariable(mesh = 1D_mesh)
>  

Re: Applying Two Different Diffusion Coefficients at Single FV Face

[Guyer, Jonathan E. Dr. (Fed)]

Yes!

> On Sep 16, 2016, at 2:18 PM, Raymond Smith  wrote:
> 
> A side note, here it may actually be more convenient to think about the D in 
> terms of the volumes, so you could also define it as a cell variable with 
> values specified as D1 when x<=1 and D2 when x>1, then use the 
> harmonicFaceValue attribute when you put it into the governing equation to 
> have FiPy do this calculation for you and avoid awkward use of the mesh 
> spacing in the specification of the variable values.
> 
> On Fri, Sep 16, 2016 at 11:12 AM, Raymond Smith  wrote:
> Hi, Ian.
> 
> I don't think there is such a thing as having two different flux coefficients 
> at the same face for the same governing PDE. The flux through a given face is 
> calculated by the coefficient at that face times some approximation of the 
> gradient in a field variable at that face, like
> D * grad(c),
> both evaluated at the face, so there is only one coefficient in that 
> expression. And in the FiPy FV approach D doesn't come into play at all 
> except at the faces, so it can't change immediately to the right or left of a 
> face but only from one face to the next.
> 
> That said, perhaps you could try following a common approach used when a FV 
> interface is directly between two bulk regions with different continuum 
> properties -- use some sort of a mean of the continuum transport coefficient 
> in the adjacent volumes as the approximation for the coefficient at the face. 
> I would suggest the harmonic mean (one simple reason: it preserves zero flux 
> through the face if either of the adjacent cells has a zero coefficient). You 
> could try:
> Dint = 2*D1*D2/(D1+D2)
> Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) & 
> (1D_mesh.faceCenters[0] < 1.0+dx/2)))
> If you do that, you should see a change in the slope of the field variable 
> right near the face between the two bulk regions, which seems to be what 
> you're missing. Of course, ideally your mesh is fine enough that the dx is 
> insignificant in terms of the solution, but this is a common approach to get 
> a slightly more reasonable estimate for what happens between two regions.
> 
> Best,
> Ray
> 
> On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian  
> wrote:
> Hello All,
> 
>  
> 
> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard elliptic 
> diffusion equation (\nabla. (D \nabla \phi) = f is being solved.
> 
>  
> 
> The diffusion coefficient, D, changes sharply at x = 1.0. Although this is 
> quite similar to examples.diffusion.mesh1D,  there is a key difference:
> 
>  
> 
> We are seeking an instantaneous change in the coefficient value at that face, 
> i.e. something of the effect that on one side of the face (i.e. from x = 0 up 
> to and including x = 1.0) diffusion coefficient D has a value D1.  Then there 
> is an abrupt change in diffusion coefficient from x >= 1.0 (very notable, 
> including the face at x = 1.0) and for the remainder of the mesh, diffusion 
> coefficient D has a value D2.  Currently, we have to either decide on whether 
> the single face at x = 1.0 gets the value D1 or D2. However, when we run with 
> either of these cases, the results are not accurate (when compared against a 
> benchmark from a commercial FV solver). We are getting a smooth interpolation 
> across the face , rather than the sharp change correctly predicted by the 
> benchmarking code.
> 
>  
> 
> With the code below, we have only been able to set the coefficient at that 
> middle face to one of either D1 or D2. This doesn’t result in the desired 
> instantaneous change in the coefficient value, but instead, the new 
> coefficient value only applies from the NEXT face (which is a distance dx 
> away) from this critical interior boundary.
> 
>  
> 
> Coefficient = FaceVariable(mesh = 1D_mesh)
>   
> # Define facevariable on mesh
> 
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0))   
>  
> # Set coefficient values in 1st half of mesh
> 
> Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) & 
> (1D_mesh.faceCenters[0] <= 2.0)))# Set coefficient values 
> in 2nd half of mesh
> 
>  
> 
> An alternative with the inequalities adjusted, but that still only permits a 
> single coefficient at that face:
> 
> Coefficient = FaceVariable(mesh = 1D_mesh)
>   
>  # Define facevariable on mesh
> 
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] < 1.0))
>   
> # Set coefficient values in 1st half of mesh
> 
> 

Re: Applying Two Different Diffusion Coefficients at Single FV Face

[Raymond Smith]

In general, if you're going to scale things in your problem, I tend to
recommend fully non-dimensionalizing a problem with reference lengths
scales, time scales, etc. and not using different references for different
regions described by the same PDE. In other words, if you have some sort of
relevant time scale, t_ref and a chosen length scale L_ref (which could be,
e.g., L1, L2, L1+L2), I would scale everything with those values, such that
D* = D/D_ref
D_ref = L_ref^2/t_ref
and so on for other quantities in the PDE's including differential
operators and field variables.
In that case, the harmonic mean of D1* and D2* is the same as the harmonic
mean of D1 and D2 divided by D_ref. If you change the references with
position along the length of your PDE, you run into issues with things like
this.

Ray

On Fri, Sep 16, 2016 at 11:36 AM, Gopalakrishnan, Krishnakumar <
krishnaku...@imperial.ac.uk> wrote:

> Hi Raymond,
>
>
>
> Thanks a lot for your inputs. We have been working along the same lines
> that you have suggested.
>
>
>
> There is a slightly more complication that wasn’t fully discussed
> (intentionally, to keep it simple, and to find out if any other options
> existed).
>
>
>
> Assuming that we are working on a normalised co-ordinate and the two bulk
> regions within the domain have non-uniform properties.  In this situation,
> we have to normalise the PDEs too.
>
>
>
> Thus, the diffusion coefficient becomes D_region/L_region  , where region
> can be either region 1 or region 2. , and L_region denotes the original
> length of the region.
>
>
>
> Now, although it’s a perfectly straightforward idea to define D in the
> cell centers and let fiPy handle the harmonic facevalue computation,  I
> don’t know if it’s meaningful to compute harmonicfacevalue of a length
> (which appears in the denominator).
>
>
>
> What do you suggest for L_region ?  Define that at the cell-centers again,
> and use arithmeticfacevalue , or stick with harmonicfacevalue ?
>
>
>
> Best Regards
>
>
>
> Krishna & Ian.
>
>
>
>
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Raymond Smith
> *Sent:* Friday, September 16, 2016 7:18 PM
> *To:* fipy@nist.gov
> *Subject:* Re: Applying Two Different Diffusion Coefficients at Single FV
> Face
>
>
>
> A side note, here it may actually be more convenient to think about the D
> in terms of the volumes, so you could also define it as a cell variable
> with values specified as D1 when x<=1 and D2 when x>1, then use the
> harmonicFaceValue attribute when you put it into the governing equation to
> have FiPy do this calculation for you and avoid awkward use of the mesh
> spacing in the specification of the variable values.
>
>
>
> On Fri, Sep 16, 2016 at 11:12 AM, Raymond Smith  wrote:
>
> Hi, Ian.
>
> I don't think there is such a thing as having two different flux
> coefficients at the same face for the same governing PDE. The flux through
> a given face is calculated by the coefficient at that face times some
> approximation of the gradient in a field variable at that face, like
> D * grad(c),
> both evaluated at the face, so there is only one coefficient in that
> expression. And in the FiPy FV approach D doesn't come into play at all
> except at the faces, so it can't change immediately to the right or left of
> a face but only from one face to the next.
>
> That said, perhaps you could try following a common approach used when a
> FV interface is directly between two bulk regions with different continuum
> properties -- use some sort of a mean of the continuum transport
> coefficient in the adjacent volumes as the approximation for the
> coefficient at the face. I would suggest the harmonic mean (one simple
> reason: it preserves zero flux through the face if either of the adjacent
> cells has a zero coefficient). You could try:
> Dint = 2*D1*D2/(D1+D2)
> Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] < 1.0+dx/2)))
> If you do that, you should see a change in the slope of the field variable
> right near the face between the two bulk regions, which seems to be what
> you're missing. Of course, ideally your mesh is fine enough that the dx is
> insignificant in terms of the solution, but this is a common approach to
> get a slightly more reasonable estimate for what happens between two
> regions.
>
>
>
> Best,
>
> Ray
>
>
>
> On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian <
> i.campbel...@imperial.ac.uk> wrote:
>
> Hello All,
>
>
>
> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard
> elliptic diffusion equation (\nabla. (D \nabla \phi) = f is being solved.
>
>
>
> The diffusion coefficient, D, changes sharply at x = 1.0. Although this is
> quite similar to examples.diffusion.mesh1D,  there is a key difference:
>
>
>
> We are seeking an instantaneous change in the coefficient value at that
> face, i.e. something of the effect that on one side of the face (i.e.
> from x = 

RE: Applying Two Different Diffusion Coefficients at Single FV Face

[Gopalakrishnan, Krishnakumar]

Hi Raymond,

Thanks a lot for your inputs. We have been working along the same lines that 
you have suggested.

There is a slightly more complication that wasn’t fully discussed 
(intentionally, to keep it simple, and to find out if any other options 
existed).

Assuming that we are working on a normalised co-ordinate and the two bulk 
regions within the domain have non-uniform properties.  In this situation, we 
have to normalise the PDEs too.

Thus, the diffusion coefficient becomes D_region/L_region  , where region can 
be either region 1 or region 2. , and L_region denotes the original length of 
the region.

Now, although it’s a perfectly straightforward idea to define D in the cell 
centers and let fiPy handle the harmonic facevalue computation,  I don’t know 
if it’s meaningful to compute harmonicfacevalue of a length (which appears in 
the denominator).

What do you suggest for L_region ?  Define that at the cell-centers again, and 
use arithmeticfacevalue , or stick with harmonicfacevalue ?

Best Regards

Krishna & Ian.



From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Raymond 
Smith
Sent: Friday, September 16, 2016 7:18 PM
To: fipy@nist.gov
Subject: Re: Applying Two Different Diffusion Coefficients at Single FV Face

A side note, here it may actually be more convenient to think about the D in 
terms of the volumes, so you could also define it as a cell variable with 
values specified as D1 when x<=1 and D2 when x>1, then use the 
harmonicFaceValue attribute when you put it into the governing equation to have 
FiPy do this calculation for you and avoid awkward use of the mesh spacing in 
the specification of the variable values.

On Fri, Sep 16, 2016 at 11:12 AM, Raymond Smith 
> wrote:
Hi, Ian.
I don't think there is such a thing as having two different flux coefficients 
at the same face for the same governing PDE. The flux through a given face is 
calculated by the coefficient at that face times some approximation of the 
gradient in a field variable at that face, like
D * grad(c),
both evaluated at the face, so there is only one coefficient in that 
expression. And in the FiPy FV approach D doesn't come into play at all except 
at the faces, so it can't change immediately to the right or left of a face but 
only from one face to the next.
That said, perhaps you could try following a common approach used when a FV 
interface is directly between two bulk regions with different continuum 
properties -- use some sort of a mean of the continuum transport coefficient in 
the adjacent volumes as the approximation for the coefficient at the face. I 
would suggest the harmonic mean (one simple reason: it preserves zero flux 
through the face if either of the adjacent cells has a zero coefficient). You 
could try:
Dint = 2*D1*D2/(D1+D2)
Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) & 
(1D_mesh.faceCenters[0] < 1.0+dx/2)))
If you do that, you should see a change in the slope of the field variable 
right near the face between the two bulk regions, which seems to be what you're 
missing. Of course, ideally your mesh is fine enough that the dx is 
insignificant in terms of the solution, but this is a common approach to get a 
slightly more reasonable estimate for what happens between two regions.

Best,
Ray

On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian 
> wrote:
Hello All,

We have a uniformly spaced 1D mesh of length 2.0, wherein a standard elliptic 
diffusion equation (\nabla. (D \nabla \phi) = f is being solved.

The diffusion coefficient, D, changes sharply at x = 1.0. Although this is 
quite similar to examples.diffusion.mesh1D,  there is a key difference:

We are seeking an instantaneous change in the coefficient value at that face, 
i.e. something of the effect that on one side of the face (i.e. from x = 0 up 
to and including x = 1.0) diffusion coefficient D has a value D1.  Then there 
is an abrupt change in diffusion coefficient from x >= 1.0 (very notable, 
including the face at x = 1.0) and for the remainder of the mesh, diffusion 
coefficient D has a value D2.  Currently, we have to either decide on whether 
the single face at x = 1.0 gets the value D1 or D2. However, when we run with 
either of these cases, the results are not accurate (when compared against a 
benchmark from a commercial FV solver). We are getting a smooth interpolation 
across the face , rather than the sharp change correctly predicted by the 
benchmarking code.

With the code below, we have only been able to set the coefficient at that 
middle face to one of either D1 or D2. This doesn’t result in the desired 
instantaneous change in the coefficient value, but instead, the new coefficient 
value only applies from the NEXT face (which is a distance dx away) from this 
critical interior boundary.

Coefficient = FaceVariable(mesh = 1D_mesh)  

Re: Applying Two Different Diffusion Coefficients at Single FV Face

[Raymond Smith]

A side note, here it may actually be more convenient to think about the D
in terms of the volumes, so you could also define it as a cell variable
with values specified as D1 when x<=1 and D2 when x>1, then use the
harmonicFaceValue attribute when you put it into the governing equation to
have FiPy do this calculation for you and avoid awkward use of the mesh
spacing in the specification of the variable values.

On Fri, Sep 16, 2016 at 11:12 AM, Raymond Smith  wrote:

> Hi, Ian.
>
> I don't think there is such a thing as having two different flux
> coefficients at the same face for the same governing PDE. The flux through
> a given face is calculated by the coefficient at that face times some
> approximation of the gradient in a field variable at that face, like
> D * grad(c),
> both evaluated at the face, so there is only one coefficient in that
> expression. And in the FiPy FV approach D doesn't come into play at all
> except at the faces, so it can't change immediately to the right or left of
> a face but only from one face to the next.
>
> That said, perhaps you could try following a common approach used when a
> FV interface is directly between two bulk regions with different continuum
> properties -- use some sort of a mean of the continuum transport
> coefficient in the adjacent volumes as the approximation for the
> coefficient at the face. I would suggest the harmonic mean (one simple
> reason: it preserves zero flux through the face if either of the adjacent
> cells has a zero coefficient). You could try:
> Dint = 2*D1*D2/(D1+D2)
> Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] < 1.0+dx/2)))
> If you do that, you should see a change in the slope of the field variable
> right near the face between the two bulk regions, which seems to be what
> you're missing. Of course, ideally your mesh is fine enough that the dx is
> insignificant in terms of the solution, but this is a common approach to
> get a slightly more reasonable estimate for what happens between two
> regions.
>
> Best,
> Ray
>
> On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian <
> i.campbel...@imperial.ac.uk> wrote:
>
>> Hello All,
>>
>>
>>
>> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard
>> elliptic diffusion equation (\nabla. (D \nabla \phi) = f is being solved.
>>
>>
>>
>> The diffusion coefficient, D, changes sharply at x = 1.0. Although this
>> is quite similar to examples.diffusion.mesh1D,  there is a key difference:
>>
>>
>>
>> We are seeking an instantaneous change in the coefficient value at that
>> face, i.e. something of the effect that on one side of the face (i.e.
>> from x = 0 up to and including x = 1.0) diffusion coefficient D has a value
>> D1.  Then there is an abrupt change in diffusion coefficient from x >= 1.0
>> (very notable, including the face at x = 1.0) and for the remainder of the
>> mesh, diffusion coefficient D has a value D2.  Currently, we have to either
>> decide on whether the single face at x = 1.0 gets the value D1 or
>> D2. However, when we run with either of these cases, the results are not
>> accurate (when compared against a benchmark from a commercial FV solver).
>> We are getting a smooth interpolation across the face , rather than the
>> sharp change correctly predicted by the benchmarking code.
>>
>>
>>
>> With the code below, we have only been able to set the coefficient at
>> that middle face to one of either D1 or D2. This doesn’t result in the
>> desired instantaneous change in the coefficient value, but instead, the new
>> coefficient value only applies from the NEXT face (which is a distance dx
>> away) from this critical interior boundary.
>>
>>
>>
>> Coefficient = FaceVariable(mesh = 1D_mesh)
>>
>>  #
>> Define facevariable on mesh
>>
>> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0))
>>
>> # Set coefficient values in 1st half of mesh
>>
>> Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) &
>> (1D_mesh.faceCenters[0] <= 2.0)))# Set coefficient
>> values in 2nd half of mesh
>>
>>
>>
>> An alternative with the inequalities adjusted, but that still only
>> permits a single coefficient at that face:
>>
>> Coefficient = FaceVariable(mesh = 1D_mesh)
>>
>>  # Define facevariable
>> on mesh
>>
>> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] < 1.0))
>>
>> # Set coefficient values in 1st half of mesh
>>
>> Coefficient.setValue(D2, where=((1.0 <= 1D_mesh.faceCenters[0]) &
>> (1D_mesh.faceCenters[0] <= 2.0)))  # Set coefficient
>> values in 2nd half of mesh
>>
>>
>>
>> Sincerely,
>>
>>
>>
>> -  Ian
>>
>>
>>
>> P.S. Daniel, thank you very much for your previous reply on sweep speeds!
>> It was very helpful.
>>
>>
>>
>> ___
>> fipy mailing list
>> 

Re: Applying Two Different Diffusion Coefficients at Single FV Face

[Raymond Smith]

Hi, Ian.

I don't think there is such a thing as having two different flux
coefficients at the same face for the same governing PDE. The flux through
a given face is calculated by the coefficient at that face times some
approximation of the gradient in a field variable at that face, like
D * grad(c),
both evaluated at the face, so there is only one coefficient in that
expression. And in the FiPy FV approach D doesn't come into play at all
except at the faces, so it can't change immediately to the right or left of
a face but only from one face to the next.

That said, perhaps you could try following a common approach used when a FV
interface is directly between two bulk regions with different continuum
properties -- use some sort of a mean of the continuum transport
coefficient in the adjacent volumes as the approximation for the
coefficient at the face. I would suggest the harmonic mean (one simple
reason: it preserves zero flux through the face if either of the adjacent
cells has a zero coefficient). You could try:
Dint = 2*D1*D2/(D1+D2)
Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) &
(1D_mesh.faceCenters[0] < 1.0+dx/2)))
If you do that, you should see a change in the slope of the field variable
right near the face between the two bulk regions, which seems to be what
you're missing. Of course, ideally your mesh is fine enough that the dx is
insignificant in terms of the solution, but this is a common approach to
get a slightly more reasonable estimate for what happens between two
regions.

Best,
Ray

On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian  wrote:

> Hello All,
>
>
>
> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard
> elliptic diffusion equation (\nabla. (D \nabla \phi) = f is being solved.
>
>
>
> The diffusion coefficient, D, changes sharply at x = 1.0. Although this is
> quite similar to examples.diffusion.mesh1D,  there is a key difference:
>
>
>
> We are seeking an instantaneous change in the coefficient value at that
> face, i.e. something of the effect that on one side of the face (i.e.
> from x = 0 up to and including x = 1.0) diffusion coefficient D has a value
> D1.  Then there is an abrupt change in diffusion coefficient from x >= 1.0
> (very notable, including the face at x = 1.0) and for the remainder of the
> mesh, diffusion coefficient D has a value D2.  Currently, we have to either
> decide on whether the single face at x = 1.0 gets the value D1 or
> D2. However, when we run with either of these cases, the results are not
> accurate (when compared against a benchmark from a commercial FV solver).
> We are getting a smooth interpolation across the face , rather than the
> sharp change correctly predicted by the benchmarking code.
>
>
>
> With the code below, we have only been able to set the coefficient at that
> middle face to one of either D1 or D2. This doesn’t result in the desired
> instantaneous change in the coefficient value, but instead, the new
> coefficient value only applies from the NEXT face (which is a distance dx
> away) from this critical interior boundary.
>
>
>
> Coefficient = FaceVariable(mesh = 1D_mesh)
>
>  #
> Define facevariable on mesh
>
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0))
>
> # Set coefficient values in 1st half of mesh
>
> Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] <= 2.0)))# Set coefficient
> values in 2nd half of mesh
>
>
>
> An alternative with the inequalities adjusted, but that still only permits
> a single coefficient at that face:
>
> Coefficient = FaceVariable(mesh = 1D_mesh)
>
>  # Define facevariable on
> mesh
>
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] < 1.0))
>
> # Set coefficient values in 1st half of mesh
>
> Coefficient.setValue(D2, where=((1.0 <= 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] <= 2.0)))  # Set coefficient
> values in 2nd half of mesh
>
>
>
> Sincerely,
>
>
>
> -  Ian
>
>
>
> P.S. Daniel, thank you very much for your previous reply on sweep speeds!
> It was very helpful.
>
>
>
> ___
> 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 ]

Applying Two Different Diffusion Coefficients at Single FV Face

[Campbell, Ian]

Hello All,

We have a uniformly spaced 1D mesh of length 2.0, wherein a standard elliptic 
diffusion equation (\nabla. (D \nabla \phi) = f is being solved.

The diffusion coefficient, D, changes sharply at x = 1.0. Although this is 
quite similar to examples.diffusion.mesh1D,  there is a key difference:

We are seeking an instantaneous change in the coefficient value at that face, 
i.e. something of the effect that on one side of the face (i.e. from x = 0 up 
to and including x = 1.0) diffusion coefficient D has a value D1.  Then there 
is an abrupt change in diffusion coefficient from x >= 1.0 (very notable, 
including the face at x = 1.0) and for the remainder of the mesh, diffusion 
coefficient D has a value D2.  Currently, we have to either decide on whether 
the single face at x = 1.0 gets the value D1 or D2. However, when we run with 
either of these cases, the results are not accurate (when compared against a 
benchmark from a commercial FV solver). We are getting a smooth interpolation 
across the face , rather than the sharp change correctly predicted by the 
benchmarking code.

With the code below, we have only been able to set the coefficient at that 
middle face to one of either D1 or D2. This doesn't result in the desired 
instantaneous change in the coefficient value, but instead, the new coefficient 
value only applies from the NEXT face (which is a distance dx away) from this 
critical interior boundary.

Coefficient = FaceVariable(mesh = 1D_mesh)  

# Define facevariable on mesh
Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0)) 
   # 
Set coefficient values in 1st half of mesh
Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) & 
(1D_mesh.faceCenters[0] <= 2.0)))# Set coefficient values 
in 2nd half of mesh

An alternative with the inequalities adjusted, but that still only permits a 
single coefficient at that face:
Coefficient = FaceVariable(mesh = 1D_mesh)  

 # Define facevariable on mesh
Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] < 1.0))  
# 
Set coefficient values in 1st half of mesh
Coefficient.setValue(D2, where=((1.0 <= 1D_mesh.faceCenters[0]) & 
(1D_mesh.faceCenters[0] <= 2.0)))  # Set coefficient values in 
2nd half of mesh

Sincerely,


-  Ian

P.S. Daniel, thank you very much for your previous reply on sweep speeds! It 
was very helpful.

___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: sources errors in advection/diffusion problems, and one solution

[James Pringle]

No worries -- I had to do it to figure out the problem in my more complex
domain and equation... The issue which surprised me was that the value the
variable was initialized to had an effect on the steady solution.

Jamie

On Fri, Sep 16, 2016 at 8:14 AM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> James -
>
> I think Daniel will have more insight into why this is happening and if
> there is anything to be done about it besides artificial relaxation, but I
> just want to say how much I appreciate your putting this together. This is
> a very lucid illustration.
>
> - Jon
>
> > On Sep 15, 2016, at 5:13 PM, James Pringle  wrote:
> >
> > Dear FiPy users --
> >
> >This is a simple example of how and why fipy may fail to solve a
> >steady advection diffusion problem, and how solving the transient
> >problem can reduce the error. I also found something that was a
> >surprise -- the "initial" condition of a steady problem can affect
> >the solution for some solvers.
> >
> >At the end are two interesting questions for those who want to
> >understand what FiPy is actually doing I admit to being a bit
> >lost
> >
> >The equation I am solving is
> >
> > \Del\dot (D\del psi + u*psi) =0
> >
> >Where the diffusion D is 1, and u is a vector (1,0) -- so there is
> >only a flow of speed -1 in the x direction.  I solve this equation
> >on a 10 by 10 grid. The boundary conditions are no normal gradient
> >on the y=0 and y=10 boundary:
> >
> > psi_y =0 at y=0 and y=10
> >
> >For the x boundary, I impose a value of x=1 on the inflow boundary at
> x=10
> >(this is a little tricky -- the way the equation is written, u is
> >the negative of velocity).
> >
> > psi=1 at x=10
> >
> >and a no-normal-gradient condition at x=0.
> >
> > psi_x=0 at x=0
> >
> >since all of the domain and boundary is symmetrical with respect to
> >y, we can assume psi_y=0 is zero everywhere. This reduces the
> equation to
> >
> > psi_xx + psi_x =0
> >
> >The general solution to this equation is
> >
> > psi=C1+C2*exp(-x)
> >
> >Where C1 and C2 are constants. For these boundary conditions, C1=1
> >and C2=0, so psi=1 everywhere.
> >
> >Now run the code SquareGrid_HomemadeDelaunay and look at figure(3)
> >-- this is the plot of psi versus x, and you can see that it does
> >not match the true solution of psi=1 everywhere! Instead, it
> >appears to be decaying exponential. The blue line is actually just
> >(1+exp(-x)). What is going on? It appears that small errors in the
> >boundary condition at x=0 are allowing C2 to not be exactly 0, and
> >this error is this exponential mode. The error is the artificial
> >exiting of a correct mode of the interior equation, albeit one that
> >should not be excited by these BC's.
> >
> >Interestingly, the size of this error depends on the value the psi
> >is initially set to. If the line
> >
> >psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
> >
> >is changed so psi is initially 1, the error goes away entirely; if
> >it is set to some other value, you get different errors. I do not
> >know if this is a bug, or just numerical error in a well programmed
> >solver.
> >
> >Now if you run SquareGrid_HomemadeDelaunay_transient  which
> implements
> >
> >  psi_t = \Del\dot (D\del psi + u*psi)
> >
> >you can see that the error in the numerical solution is advected
> >out of the x=0 boundary, and the solution approaches the true
> >solution of psi=1 rapidly.
> >
> >The interesting question is if the formulation of the boundary
> >condition at x=0 could be altered to less excite the spurious mode?
> >
> >Also, why does the "initial" condition have any effect on the
> >steady equation?
> >
> >Cheers,
> >Jamie
> >
> >  transient.py>__
> _
> > fipy mailing list
> > fipy@nist.gov
> > https://urldefense.proofpoint.com/v2/url?u=http-3A__www.
> ctcms.nist.gov_fipy=DQICAg=c6MrceVCY5m5A_KAUkrdoA=
> 7HJI3EH6RqKWf16fbYIYKw=IaTL4nEetwjm4G4qfj8cwLzvztKzui
> Fsw_Nhksv_oWQ=vG-nxTf76KxE_CqEHTjt2jIkoy0l9M6X8bm01ypXaBQ=
> >  [ NIST internal ONLY: https://urldefense.proofpoint.
> com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_
> fipy=DQICAg=c6MrceVCY5m5A_KAUkrdoA=7HJI3EH6RqKWf16fbYIYKw=
> IaTL4nEetwjm4G4qfj8cwLzvztKzuiFsw_Nhksv_oWQ=
> 7JGh0Zz82O69cJiLFKmkV3NfW2TLz6KB_ngkAGCrYGI=  ]
>
>
> ___
> fipy mailing list
> fipy@nist.gov
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.
> ctcms.nist.gov_fipy=DQICAg=c6MrceVCY5m5A_KAUkrdoA=
> 7HJI3EH6RqKWf16fbYIYKw=IaTL4nEetwjm4G4qfj8cwLzvztKzui
> Fsw_Nhksv_oWQ=vG-nxTf76KxE_CqEHTjt2jIkoy0l9M6X8bm01ypXaBQ=
>   [ NIST internal ONLY: https://urldefense.proofpoint.
> com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_

Re: sources errors in advection/diffusion problems, and one solution

[Guyer, Jonathan E. Dr. (Fed)]

James -

I think Daniel will have more insight into why this is happening and if there 
is anything to be done about it besides artificial relaxation, but I just want 
to say how much I appreciate your putting this together. This is a very lucid 
illustration.

- Jon

> On Sep 15, 2016, at 5:13 PM, James Pringle  wrote:
> 
> Dear FiPy users --
> 
>This is a simple example of how and why fipy may fail to solve a
>steady advection diffusion problem, and how solving the transient
>problem can reduce the error. I also found something that was a
>surprise -- the "initial" condition of a steady problem can affect
>the solution for some solvers.
> 
>At the end are two interesting questions for those who want to
>understand what FiPy is actually doing I admit to being a bit
>lost
> 
>The equation I am solving is
> 
> \Del\dot (D\del psi + u*psi) =0
> 
>Where the diffusion D is 1, and u is a vector (1,0) -- so there is
>only a flow of speed -1 in the x direction.  I solve this equation
>on a 10 by 10 grid. The boundary conditions are no normal gradient
>on the y=0 and y=10 boundary:
> 
> psi_y =0 at y=0 and y=10
> 
>For the x boundary, I impose a value of x=1 on the inflow boundary at x=10
>(this is a little tricky -- the way the equation is written, u is
>the negative of velocity).
> 
> psi=1 at x=10
> 
>and a no-normal-gradient condition at x=0.
> 
> psi_x=0 at x=0
> 
>since all of the domain and boundary is symmetrical with respect to
>y, we can assume psi_y=0 is zero everywhere. This reduces the equation to
> 
> psi_xx + psi_x =0
> 
>The general solution to this equation is
> 
> psi=C1+C2*exp(-x)
> 
>Where C1 and C2 are constants. For these boundary conditions, C1=1
>and C2=0, so psi=1 everywhere.
> 
>Now run the code SquareGrid_HomemadeDelaunay and look at figure(3)
>-- this is the plot of psi versus x, and you can see that it does
>not match the true solution of psi=1 everywhere! Instead, it
>appears to be decaying exponential. The blue line is actually just
>(1+exp(-x)). What is going on? It appears that small errors in the
>boundary condition at x=0 are allowing C2 to not be exactly 0, and
>this error is this exponential mode. The error is the artificial
>exiting of a correct mode of the interior equation, albeit one that
>should not be excited by these BC's.
> 
>Interestingly, the size of this error depends on the value the psi
>is initially set to. If the line
> 
>psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
> 
>is changed so psi is initially 1, the error goes away entirely; if
>it is set to some other value, you get different errors. I do not
>know if this is a bug, or just numerical error in a well programmed
>solver.
> 
>Now if you run SquareGrid_HomemadeDelaunay_transient  which implements
> 
>  psi_t = \Del\dot (D\del psi + u*psi) 
> 
>you can see that the error in the numerical solution is advected
>out of the x=0 boundary, and the solution approaches the true
>solution of psi=1 rapidly.
> 
>The interesting question is if the formulation of the boundary
>condition at x=0 could be altered to less excite the spurious mode?
> 
>Also, why does the "initial" condition have any effect on the
>steady equation?
> 
>Cheers,
>Jamie
> 
> ___
> 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 ]