Re: IndexError: index is out of bounds for axis 1 with size

[Daniel Wheeler]

One possible way to handle this would be to have Pu and Pd defined
everywhere, but only have meaning in the edge cells. The coefficients C3
and C4 would only have a non-zero value in the end cells. Move the end
cells so the cell centers are on the boundaries. The boundary conditions
for the Pf and Pk equations can then be source terms in the edge cells that
constrain the edge values. The equations can then be tightly coupled
implicitly. I think this will all work quite well. The important thing to
know is that a single cell value can be constrained with the application of
a large implicit and explicit source, e.g.,

x, y = mesh.cellCenters
big_coeff = CellVariable(mesh=mesh, value=0.)
big_coeff[x < dx] = 1e+20
big_coeff[x > L - dx] = 1e+20
TransientTerm() == big_coeff * value - ImplicitSourceTerm(big_coeff)

if `L` is the length of the mesh and `dx` is the grid spacing. The edge
cells will be constrained to `value`. `value` can be one of the
CellVariables that are in the other equations.

I hope that helps you get started.

On Mon, May 9, 2016 at 11:57 AM, Mohammad Kazemi  wrote:

> I forgot to mention that the [image: P_f] and [image: P_k] values in eqpd
> and eqpu are actually as follow,
>
> eqpd:
>
> [image: -\frac{\partial P_d}{\partial t} = C_3 (P_f (0,t) - P_d(t)) + C_4
> (P_k(0,t)-P_d(t))]
>
> and eqpu:
>
> [image: \frac{\partial P_u}{\partial t} = C_3 (P_f (L,t) - P_u(t)) + C_4
> (P_k(L,t)-P_u(t))]
>
>
> On Mon, May 9, 2016 at 11:19 AM, Mohammad Kazemi 
> wrote:
>
>> Thanks for your response. Here are my equations:
>>
>> eqpk:
>>
>> [image: \frac{\partial P_k}{\partial t} = D \frac{\partial^2
>> P_k}{\partial x^2}]
>> with IC and BCs:
>> [image: P_k (x,0) = Pi]
>> [image: P_k(0,t) = P_d(t)]
>> [image: P_k(L,t) = P_u(t)]
>>
>> eqpf:
>>
>> ∂Pf/∂t = C1 ∂2 Pf/∂x2 ​ - C2 (Pf-Pk)
>>
>> with IC and BCs:
>>
>> [image: P_f (x,0) = Pi]
>> [image: P_f(0,t) = P_d(t)]
>> [image: P_f(L,t) = P_u(t)]
>>
>> eqpd:
>>
>> [image: - \frac{\partial P_d}{\partial t} = C_3 (P_f - P_k)+C_4 (P_k-P_d)]
>>
>> with IC:
>>
>> [image: P_d (0) = Pi]
>>
>> eqpu:
>>
>> [image: \frac{\partial P_u}{\partial t} = C_3 (P_f - P_u)+C_4 (P_k-P_u)]
>>
>> with IC:
>>
>> [image: P_u (0) = Pi]
>>
>> So [image: P_k] and [image: P_f] are functions of space and time while 
>> [image:
>> P_d] and [image: P_u] are only function of time. How should i define a
>> variable which is only a function of time? If I want to have schematic of
>> problem, it is as below.
>>
>>
>>
>> ​
>>
>> The pressure in tanks [image: P_u] and [image: P_d] are uniform and only
>> changes with time.
>>
>> I appreciate it,
>>
>>
>>
>>
>> On Mon, May 9, 2016 at 9:43 AM, Guyer, Jonathan E. Dr. (Fed) <
>> jonathan.gu...@nist.gov> wrote:
>>
>>> By declaring Pd and Pu as CellVariables, you are declaring them to be
>>> functions of space. That's what a CellVariable is.
>>>
>>> Furthermore, eqpd and eqpu depend on Pf and Pk, so that necessitates
>>> that Pd and Pu are functions of space.
>>>
>>> Can you show us the math for the equations you are trying to solve?
>>>
>>> > On May 7, 2016, at 1:12 PM, Mohammad Kazemi 
>>> wrote:
>>> >
>>> > Thank you. It is running now, but I have a question. Pf and Pk are
>>> functions of location (x) and time (t). However, Pd and Pu are only
>>> functions of time (they represent pressure at downstream and upstream
>>> tanks). When I solve this problem, I will get a vector for Pd and Pu (at
>>> last timestep) which varies with location. For example, for Pu,
>>> >
>>> > print Pu
>>> >
>>> > [ 7171802.05029459  7171719.44929532  7171632.68237858
>>> 7171538.95796536
>>> >   7171435.46595385  7171319.37283423  7171187.81897787
>>> 7171037.91833804
>>> >   7170866.76056175  7170671.41527712]
>>> >
>>> > Which shows the Pu values at 10 different mesh blocks. Is it related
>>> to how I define my cell variables?
>>> >
>>> > I appreciate your helps.
>>> >
>>> > Bests,
>>> >
>>> > On Fri, May 6, 2016 at 5:21 PM, Daniel Wheeler <
>>> daniel.wheel...@gmail.com> wrote:
>>> > The following runs for me with a few changes.
>>> >
>>> > On Fri, May 6, 2016 at 1:25 PM, Mohammad Kazemi 
>>> wrote:
>>> > >
>>> > > Pk.constrain([Pd], mesh.facesLeft)
>>> > > Pk.constrain([Pu], mesh.facesRight)
>>> >
>>> > Pk.constrain(Pd.faceValue, mesh.facesLeft)
>>> > Pk.constrain(Pu.faceValue, mesh.facesRight)
>>> >
>>> > > Pf.constrain([Pd], mesh.facesLeft)
>>> > > Pf.constrain([Pu], mesh.facesRight)
>>> >
>>> > Pf.constrain(Pd.faceValue, mesh.facesLeft)
>>> > Pf.constrain(Pu.faceValue, mesh.facesRight)
>>> >
>>> > This makes sense since the Pf and Pk are scalars and [Pd] and [Pu] are
>>> > effectively vectors. Also the face value is required since the
>>> > constraint is on the faces.
>>> >
>>> > I hope this helps.
>>> >
>>> > --
>>> > 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: IndexError: index is out of bounds for axis 1 with size

[Mohammad Kazemi]

I forgot to mention that the [image: P_f] and [image: P_k] values in eqpd
and eqpu are actually as follow,

eqpd:

[image: -\frac{\partial P_d}{\partial t} = C_3 (P_f (0,t) - P_d(t)) + C_4
(P_k(0,t)-P_d(t))]

and eqpu:

[image: \frac{\partial P_u}{\partial t} = C_3 (P_f (L,t) - P_u(t)) + C_4
(P_k(L,t)-P_u(t))]


On Mon, May 9, 2016 at 11:19 AM, Mohammad Kazemi  wrote:

> Thanks for your response. Here are my equations:
>
> eqpk:
>
> [image: \frac{\partial P_k}{\partial t} = D \frac{\partial^2 P_k}{\partial
> x^2}]
> with IC and BCs:
> [image: P_k (x,0) = Pi]
> [image: P_k(0,t) = P_d(t)]
> [image: P_k(L,t) = P_u(t)]
>
> eqpf:
>
> ∂Pf/∂t = C1 ∂2 Pf/∂x2 ​ - C2 (Pf-Pk)
>
> with IC and BCs:
>
> [image: P_f (x,0) = Pi]
> [image: P_f(0,t) = P_d(t)]
> [image: P_f(L,t) = P_u(t)]
>
> eqpd:
>
> [image: - \frac{\partial P_d}{\partial t} = C_3 (P_f - P_k)+C_4 (P_k-P_d)]
>
> with IC:
>
> [image: P_d (0) = Pi]
>
> eqpu:
>
> [image: \frac{\partial P_u}{\partial t} = C_3 (P_f - P_u)+C_4 (P_k-P_u)]
>
> with IC:
>
> [image: P_u (0) = Pi]
>
> So [image: P_k] and [image: P_f] are functions of space and time while [image:
> P_d] and [image: P_u] are only function of time. How should i define a
> variable which is only a function of time? If I want to have schematic of
> problem, it is as below.
>
>
>
> ​
>
> The pressure in tanks [image: P_u] and [image: P_d] are uniform and only
> changes with time.
>
> I appreciate it,
>
>
>
>
> On Mon, May 9, 2016 at 9:43 AM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
>
>> By declaring Pd and Pu as CellVariables, you are declaring them to be
>> functions of space. That's what a CellVariable is.
>>
>> Furthermore, eqpd and eqpu depend on Pf and Pk, so that necessitates that
>> Pd and Pu are functions of space.
>>
>> Can you show us the math for the equations you are trying to solve?
>>
>> > On May 7, 2016, at 1:12 PM, Mohammad Kazemi  wrote:
>> >
>> > Thank you. It is running now, but I have a question. Pf and Pk are
>> functions of location (x) and time (t). However, Pd and Pu are only
>> functions of time (they represent pressure at downstream and upstream
>> tanks). When I solve this problem, I will get a vector for Pd and Pu (at
>> last timestep) which varies with location. For example, for Pu,
>> >
>> > print Pu
>> >
>> > [ 7171802.05029459  7171719.44929532  7171632.68237858  7171538.95796536
>> >   7171435.46595385  7171319.37283423  7171187.81897787  7171037.91833804
>> >   7170866.76056175  7170671.41527712]
>> >
>> > Which shows the Pu values at 10 different mesh blocks. Is it related to
>> how I define my cell variables?
>> >
>> > I appreciate your helps.
>> >
>> > Bests,
>> >
>> > On Fri, May 6, 2016 at 5:21 PM, Daniel Wheeler <
>> daniel.wheel...@gmail.com> wrote:
>> > The following runs for me with a few changes.
>> >
>> > On Fri, May 6, 2016 at 1:25 PM, Mohammad Kazemi 
>> wrote:
>> > >
>> > > Pk.constrain([Pd], mesh.facesLeft)
>> > > Pk.constrain([Pu], mesh.facesRight)
>> >
>> > Pk.constrain(Pd.faceValue, mesh.facesLeft)
>> > Pk.constrain(Pu.faceValue, mesh.facesRight)
>> >
>> > > Pf.constrain([Pd], mesh.facesLeft)
>> > > Pf.constrain([Pu], mesh.facesRight)
>> >
>> > Pf.constrain(Pd.faceValue, mesh.facesLeft)
>> > Pf.constrain(Pu.faceValue, mesh.facesRight)
>> >
>> > This makes sense since the Pf and Pk are scalars and [Pd] and [Pu] are
>> > effectively vectors. Also the face value is required since the
>> > constraint is on the faces.
>> >
>> > I hope this helps.
>> >
>> > --
>> > Daniel Wheeler
>> > ___
>> > fipy mailing list
>> > fipy@nist.gov
>> > http://www.ctcms.nist.gov/fipy
>> >   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>> >
>> >
>> >
>> > --
>> > Mohammad Kazemi
>> > West Virginia University
>> > Department of Petroleum and Natural Gas Engineering
>> > ___
>> > 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 ]
>>
>
>
>
> --
> Mohammad Kazemi
> West Virginia University
> Department of Petroleum and Natural Gas Engineering
>



-- 
Mohammad Kazemi
West Virginia University
Department of Petroleum and Natural Gas Engineering
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: IndexError: index is out of bounds for axis 1 with size

[Mohammad Kazemi]

Thanks for your response. Here are my equations:

eqpk:

[image: \frac{\partial P_k}{\partial t} = D \frac{\partial^2 P_k}{\partial
x^2}]
with IC and BCs:
[image: P_k (x,0) = Pi]
[image: P_k(0,t) = P_d(t)]
[image: P_k(L,t) = P_u(t)]

eqpf:

∂Pf/∂t = C1 ∂2 Pf/∂x2 ​ - C2 (Pf-Pk)

with IC and BCs:

[image: P_f (x,0) = Pi]
[image: P_f(0,t) = P_d(t)]
[image: P_f(L,t) = P_u(t)]

eqpd:

[image: - \frac{\partial P_d}{\partial t} = C_3 (P_f - P_k)+C_4 (P_k-P_d)]

with IC:

[image: P_d (0) = Pi]

eqpu:

[image: \frac{\partial P_u}{\partial t} = C_3 (P_f - P_u)+C_4 (P_k-P_u)]

with IC:

[image: P_u (0) = Pi]

So [image: P_k] and [image: P_f] are functions of space and time while [image:
P_d] and [image: P_u] are only function of time. How should i define a
variable which is only a function of time? If I want to have schematic of
problem, it is as below.



​

The pressure in tanks [image: P_u] and [image: P_d] are uniform and only
changes with time.

I appreciate it,




On Mon, May 9, 2016 at 9:43 AM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> By declaring Pd and Pu as CellVariables, you are declaring them to be
> functions of space. That's what a CellVariable is.
>
> Furthermore, eqpd and eqpu depend on Pf and Pk, so that necessitates that
> Pd and Pu are functions of space.
>
> Can you show us the math for the equations you are trying to solve?
>
> > On May 7, 2016, at 1:12 PM, Mohammad Kazemi  wrote:
> >
> > Thank you. It is running now, but I have a question. Pf and Pk are
> functions of location (x) and time (t). However, Pd and Pu are only
> functions of time (they represent pressure at downstream and upstream
> tanks). When I solve this problem, I will get a vector for Pd and Pu (at
> last timestep) which varies with location. For example, for Pu,
> >
> > print Pu
> >
> > [ 7171802.05029459  7171719.44929532  7171632.68237858  7171538.95796536
> >   7171435.46595385  7171319.37283423  7171187.81897787  7171037.91833804
> >   7170866.76056175  7170671.41527712]
> >
> > Which shows the Pu values at 10 different mesh blocks. Is it related to
> how I define my cell variables?
> >
> > I appreciate your helps.
> >
> > Bests,
> >
> > On Fri, May 6, 2016 at 5:21 PM, Daniel Wheeler <
> daniel.wheel...@gmail.com> wrote:
> > The following runs for me with a few changes.
> >
> > On Fri, May 6, 2016 at 1:25 PM, Mohammad Kazemi 
> wrote:
> > >
> > > Pk.constrain([Pd], mesh.facesLeft)
> > > Pk.constrain([Pu], mesh.facesRight)
> >
> > Pk.constrain(Pd.faceValue, mesh.facesLeft)
> > Pk.constrain(Pu.faceValue, mesh.facesRight)
> >
> > > Pf.constrain([Pd], mesh.facesLeft)
> > > Pf.constrain([Pu], mesh.facesRight)
> >
> > Pf.constrain(Pd.faceValue, mesh.facesLeft)
> > Pf.constrain(Pu.faceValue, mesh.facesRight)
> >
> > This makes sense since the Pf and Pk are scalars and [Pd] and [Pu] are
> > effectively vectors. Also the face value is required since the
> > constraint is on the faces.
> >
> > I hope this helps.
> >
> > --
> > Daniel Wheeler
> > ___
> > fipy mailing list
> > fipy@nist.gov
> > http://www.ctcms.nist.gov/fipy
> >   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> >
> >
> >
> > --
> > Mohammad Kazemi
> > West Virginia University
> > Department of Petroleum and Natural Gas Engineering
> > ___
> > 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 ]
>



-- 
Mohammad Kazemi
West Virginia University
Department of Petroleum and Natural Gas Engineering
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: IndexError: index is out of bounds for axis 1 with size

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

By declaring Pd and Pu as CellVariables, you are declaring them to be functions 
of space. That's what a CellVariable is.

Furthermore, eqpd and eqpu depend on Pf and Pk, so that necessitates that Pd 
and Pu are functions of space. 

Can you show us the math for the equations you are trying to solve?

> On May 7, 2016, at 1:12 PM, Mohammad Kazemi  wrote:
> 
> Thank you. It is running now, but I have a question. Pf and Pk are functions 
> of location (x) and time (t). However, Pd and Pu are only functions of time 
> (they represent pressure at downstream and upstream tanks). When I solve this 
> problem, I will get a vector for Pd and Pu (at last timestep) which varies 
> with location. For example, for Pu,
> 
> print Pu
> 
> [ 7171802.05029459  7171719.44929532  7171632.68237858  7171538.95796536
>   7171435.46595385  7171319.37283423  7171187.81897787  7171037.91833804
>   7170866.76056175  7170671.41527712]
> 
> Which shows the Pu values at 10 different mesh blocks. Is it related to how I 
> define my cell variables?
> 
> I appreciate your helps.
> 
> Bests,
> 
> On Fri, May 6, 2016 at 5:21 PM, Daniel Wheeler  
> wrote:
> The following runs for me with a few changes.
> 
> On Fri, May 6, 2016 at 1:25 PM, Mohammad Kazemi  wrote:
> >
> > Pk.constrain([Pd], mesh.facesLeft)
> > Pk.constrain([Pu], mesh.facesRight)
> 
> Pk.constrain(Pd.faceValue, mesh.facesLeft)
> Pk.constrain(Pu.faceValue, mesh.facesRight)
> 
> > Pf.constrain([Pd], mesh.facesLeft)
> > Pf.constrain([Pu], mesh.facesRight)
> 
> Pf.constrain(Pd.faceValue, mesh.facesLeft)
> Pf.constrain(Pu.faceValue, mesh.facesRight)
> 
> This makes sense since the Pf and Pk are scalars and [Pd] and [Pu] are
> effectively vectors. Also the face value is required since the
> constraint is on the faces.
> 
> I hope this helps.
> 
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> 
> 
> 
> -- 
> Mohammad Kazemi
> West Virginia University
> Department of Petroleum and Natural Gas Engineering
> ___
> 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: IndexError: index is out of bounds for axis 1 with size

[Mohammad Kazemi]

Thank you. It is running now, but I have a question. Pf and Pk are
functions of location (x) and time (t). However, Pd and Pu are only
functions of time (they represent pressure at downstream and upstream
tanks). When I solve this problem, I will get a vector for Pd and Pu (at
last timestep) which varies with location. For example, for Pu,

print Pu

[ 7171802.05029459  7171719.44929532  7171632.68237858  7171538.95796536
  7171435.46595385  7171319.37283423  7171187.81897787  7171037.91833804
  7170866.76056175  7170671.41527712]

Which shows the Pu values at 10 different mesh blocks. Is it related to how
I define my cell variables?

I appreciate your helps.

Bests,

On Fri, May 6, 2016 at 5:21 PM, Daniel Wheeler 
wrote:

> The following runs for me with a few changes.
>
> On Fri, May 6, 2016 at 1:25 PM, Mohammad Kazemi 
> wrote:
> >
> > Pk.constrain([Pd], mesh.facesLeft)
> > Pk.constrain([Pu], mesh.facesRight)
>
> Pk.constrain(Pd.faceValue, mesh.facesLeft)
> Pk.constrain(Pu.faceValue, mesh.facesRight)
>
> > Pf.constrain([Pd], mesh.facesLeft)
> > Pf.constrain([Pu], mesh.facesRight)
>
> Pf.constrain(Pd.faceValue, mesh.facesLeft)
> Pf.constrain(Pu.faceValue, mesh.facesRight)
>
> This makes sense since the Pf and Pk are scalars and [Pd] and [Pu] are
> effectively vectors. Also the face value is required since the
> constraint is on the faces.
>
> I hope this helps.
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>



-- 
Mohammad Kazemi
West Virginia University
Department of Petroleum and Natural Gas Engineering
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: IndexError: index is out of bounds for axis 1 with size

[Daniel Wheeler]

The following runs for me with a few changes.

On Fri, May 6, 2016 at 1:25 PM, Mohammad Kazemi  wrote:
>
> Pk.constrain([Pd], mesh.facesLeft)
> Pk.constrain([Pu], mesh.facesRight)

Pk.constrain(Pd.faceValue, mesh.facesLeft)
Pk.constrain(Pu.faceValue, mesh.facesRight)

> Pf.constrain([Pd], mesh.facesLeft)
> Pf.constrain([Pu], mesh.facesRight)

Pf.constrain(Pd.faceValue, mesh.facesLeft)
Pf.constrain(Pu.faceValue, mesh.facesRight)

This makes sense since the Pf and Pk are scalars and [Pd] and [Pu] are
effectively vectors. Also the face value is required since the
constraint is on the faces.

I hope this helps.

-- 
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: IndexError: index is out of bounds for axis 1 with size

[Mohammad Kazemi]

Sorry for the confusion. For some reasons, it did not show me that error,
but I fixed the Pf definition. Here is my code which still give me the
index error.

CODE

from __future__ import division
from numpy import *
from pylab import *
from matplotlib import rc
from fipy import *

 INITIALLIZATION
###
L = 0.0508  #+0.0021+0.0014
nx = 10.0
dx = L / nx
Pdi=6894757.29
Pui=7170547.58
D=1e-6
tmax = 1000
sigma = 400
k0 = 10e-16
T = 300
M = 0.0399
R = 8.314
muf=1.96e-5
hf=30e-6
phif=0.25
A1=1.14
A2=0.28
wkf=4
wkd=4
wku=4

mesh = Grid1D(nx=nx, dx=dx)

### VARIABLE DEFINITIN
###
Pk = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pf = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pd = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pu = CellVariable(mesh=mesh, value=Pui,hasOld=True)
Kn=sqrt((3.14*R*T/(2*M))*muf/(Pf*hf))
# Kerogen mass balance

Pk.constrain([Pd], mesh.facesLeft)
Pk.constrain([Pu], mesh.facesRight)
eqpk = TransientTerm(var=Pk) == DiffusionTerm(D, var=Pk)

# Natural Fracture mass balance
Df = k0 * (1 + 6 * A1 * Kn + 12 * A2 * Kn**2)*Pf/muf

wkfD=wkf*D
Pf.constrain([Pd], mesh.facesLeft)
Pf.constrain([Pu], mesh.facesRight)
eqpf = TransientTerm(phif,var=Pf) == DiffusionTerm(Df,var=Pf) -
ImplicitSourceTerm(coeff=wkfD,var=Pf) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)

 Downstream tank
sigmaD = (k0*(1 + 6 * A1 * Kn + 12 * A2 * Kn**2)*sigma)/(muf)

eqpd = TransientTerm(var=Pd) == ImplicitSourceTerm(coeff=sigmaD,var=Pf) -
ImplicitSourceTerm(coeff=sigmaD,var=Pd)-
ImplicitSourceTerm(coeff=wkfD,var=Pd) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)

##  Upstream tank

eqpu = - TransientTerm(var=Pu) == ImplicitSourceTerm(coeff=sigmaD,var=Pf) -
ImplicitSourceTerm(coeff=sigmaD,var=Pu)-
ImplicitSourceTerm(coeff=wkfD,var=Pu) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)

## TIMESTEP
dt = 0.9 * dx**2 / (2 * D)
steps = 100

# Coupling

eqn = eqpk & eqpf & eqpd & eqpu

 Solver
for t in range(1):
Pk.updateOld()
Pf.updateOld()
Pd.updateOld()
Pu.updateOld()
eqn.sweep(dt=dt)



--
Error:

--

--IndexError
   Traceback (most recent call
last) in () 72
Pd.updateOld() 73 Pu.updateOld()---> 74 eqn.sweep(dt=dt)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/term.pyc
in sweep(self, var, solver, boundaryConditions, dt, underRelaxation,
residualFn, cacheResidual, cacheError)234 235 """-->
236 solver = self._prepareLinearSystem(var=var, solver=solver,
boundaryConditions=boundaryConditions, dt=dt)237
solver._applyUnderRelaxation(underRelaxation=underRelaxation)238
  residual = solver._calcResidual(residualFn=residualFn)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/term.pyc
in _prepareLinearSystem(self, var, solver, boundaryConditions, dt)
168
transientGeomCoeff=self._getTransientGeomCoeff(var),169

diffusionGeomCoeff=self._getDiffusionGeomCoeff(var),--> 170

buildExplicitIfOther=self._buildExplcitIfOther)171 172
self._buildCache(matrix, RHSvector)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/coupledBinaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
120

transientGeomCoeff=uncoupledTerm._getTransientGeomCoeff(tmpVar),
121

diffusionGeomCoeff=uncoupledTerm._getDiffusionGeomCoeff(tmpVar),-->
122
  buildExplicitIfOther=buildExplicitIfOther)123
 124 termMatrix += tmpMatrix
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/binaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
66
transientGeomCoeff=transientGeomCoeff, 67

diffusionGeomCoeff=diffusionGeomCoeff,---> 68

buildExplicitIfOther=buildExplicitIfOther) 69  70
matrix += tmpMatrix
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/unaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
97dt=dt,
98
transientGeomCoeff=transientGeomCoeff,---> 99
  diffusionGeomCoeff=diffusionGeomCoeff)
 100 elif buildExplicitIfOther:101 _, matrix,
RHSvector = self._buildMatrix(self.var,
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/abstractDiffusionTerm.pyc
in _buildMatrix(s

Re: IndexError: index is out of bounds for axis 1 with size

[Raymond Smith]

I think Daniel's point is that Pf is not defined yet when it is first
called in the script. You could fix this by moving the line using Pf that
Daniel pointed to after the definition that you refer to. I think the
general point, however, is that it makes it much easier to help if you try
to make sure that the script you send gives the errors that you see.

Best,
Ray


Virus-free.
www.avast.com

<#m_-68387843283325897_DDB4FAA8-2DD7-40BB-A1B8-4E2AA1F9FDF2>

On Fri, May 6, 2016 at 1:00 PM, Mohammad Kazemi  wrote:

> Thanks Daniel for your response. I thought defining variables as
> cellvariable should be enough.
>
> Pf = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
>
>
> Should I define it in some other way?
>
> Thanks
>
> On Fri, May 6, 2016 at 12:10 PM, Daniel Wheeler  > wrote:
>
>> On Fri, May 6, 2016 at 12:09 PM, Daniel Wheeler
>>  wrote:
>> > Hi Mo,
>>
>> Apologies for that.
>>
>> --
>> Daniel Wheeler
>> ___
>> fipy mailing list
>> fipy@nist.gov
>> http://www.ctcms.nist.gov/fipy
>>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>
>
>
>
> --
> Mohammad Kazemi
> West Virginia University
> Department of Petroleum and Natural Gas Engineering
>
> ___
> 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: IndexError: index is out of bounds for axis 1 with size

[Mohammad Kazemi]

Thanks Daniel for your response. I thought defining variables as
cellvariable should be enough.

Pf = CellVariable(mesh=mesh, value=Pdi,hasOld=True)


Should I define it in some other way?

Thanks

On Fri, May 6, 2016 at 12:10 PM, Daniel Wheeler 
wrote:

> On Fri, May 6, 2016 at 12:09 PM, Daniel Wheeler
>  wrote:
> > Hi Mo,
>
> Apologies for that.
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>



-- 
Mohammad Kazemi
West Virginia University
Department of Petroleum and Natural Gas Engineering
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: IndexError: index is out of bounds for axis 1 with size

[Daniel Wheeler]

On Fri, May 6, 2016 at 12:09 PM, Daniel Wheeler
 wrote:
> Hi Mo,

Apologies for that.

-- 
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: IndexError: index is out of bounds for axis 1 with size

[Daniel Wheeler]

Hi Mo,

I tried running your code but I got "name 'Pf' is not defined" -- a
different error.

On Fri, May 6, 2016 at 12:41 AM, Mohammad Kazemi  wrote:
>
> CODE:
> 
> from __future__ import division
> from numpy import *
> from pylab import *
> from matplotlib import rc
> from fipy import *
>
>  INITIALLIZATION
> ###
> L = 0.0508  #+0.0021+0.0014
> nx = 10.0
> dx = L / nx
> Pdi=6894757.29
> Pui=7170547.58
> D=1e-6
> tmax = 1000
> sigma = 400
> k0 = 10e-16
> T = 300
> M = 0.0399
> R = 8.314
> muf=1.96e-5
> hf=30e-6
> phif=0.25
> A1=1.14
> A2=0.28
> wkf=4
> wkd=4
> wku=4
> Kn=sqrt((3.14*R*T/(2*M))*muf/(Pf*hf))

At this point "Pf" is not defined.

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

IndexError: index is out of bounds for axis 1 with size

[Mohammad Kazemi]

Dear Fipy users and developers,

I was trying to solve a system of partial differential equations using fipy
but I am keep getting error of "" . I would appreciate if you could help me
on this error. Here is my code and error.


CODE:

from __future__ import division
from numpy import *
from pylab import *
from matplotlib import rc
from fipy import *

 INITIALLIZATION
###
L = 0.0508  #+0.0021+0.0014
nx = 10.0
dx = L / nx
Pdi=6894757.29
Pui=7170547.58
D=1e-6
tmax = 1000
sigma = 400
k0 = 10e-16
T = 300
M = 0.0399
R = 8.314
muf=1.96e-5
hf=30e-6
phif=0.25
A1=1.14
A2=0.28
wkf=4
wkd=4
wku=4
Kn=sqrt((3.14*R*T/(2*M))*muf/(Pf*hf))
mesh = Grid1D(nx=nx, dx=dx)
X = mesh.faceCenters[0]
#meshd = Grid1D(nx=nx, dx=dx)
### VARIABLE DEFINITION
Pk = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pf = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pd = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pu = CellVariable(mesh=mesh, value=Pui,hasOld=True)

# Kerogen mass balance

Pk.constrain([Pd], mesh.facesLeft)
Pk.constrain([Pu], mesh.facesRight)
eqpk = TransientTerm(var=Pk) == DiffusionTerm(D, var=Pk)

# Natural Fracture mass balance
Df = k0 * (1 + 6 * A1 * Kn + 12 * A2 * Kn**2)*Pf/muf

wkfD=wkf*D
Pf.constrain([Pd], mesh.facesLeft)
Pf.constrain([Pu], mesh.facesRight)
eqpf = TransientTerm(phif,var=Pf) == DiffusionTerm(Df,var=Pf) -
ImplicitSourceTerm(coeff=wkfD,var=Pf) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)

 Downstream tank
sigmaD = (k0*(1 + 6 * A1 * Kn + 12 * A2 * Kn**2)*sigma)/(muf)

eqpd = TransientTerm(var=Pd) == ImplicitSourceTerm(coeff=sigmaD,var=Pf) -
ImplicitSourceTerm(coeff=sigmaD,var=Pd)-
ImplicitSourceTerm(coeff=wkfD,var=Pd) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)

##  Upstream tank

eqpu = - TransientTerm(var=Pu) == ImplicitSourceTerm(coeff=sigmaD,var=Pf) -
ImplicitSourceTerm(coeff=sigmaD,var=Pu)-
ImplicitSourceTerm(coeff=wkfD,var=Pu) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)

## TIMESTEP
dt = 0.9 * dx**2 / (2 * D)
steps = 100

# Coupling

eqn = eqpk & eqpf & eqpd & eqpu

 Solver
for t in range(1):
Pk.updateOld()
Pf.updateOld()
Pd.updateOld()
Pu.updateOld()
eqn.sweep(dt=dt)

-

ERROR


IndexErrorTraceback (most recent call
last) in () 73
Pd.updateOld() 74 Pu.updateOld()---> 75 eqn.sweep(dt=dt)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/term.pyc
in sweep(self, var, solver, boundaryConditions, dt, underRelaxation,
residualFn, cacheResidual, cacheError)234 235 """-->
236 solver = self._prepareLinearSystem(var=var, solver=solver,
boundaryConditions=boundaryConditions, dt=dt)237
solver._applyUnderRelaxation(underRelaxation=underRelaxation)238
  residual = solver._calcResidual(residualFn=residualFn)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/term.pyc
in _prepareLinearSystem(self, var, solver, boundaryConditions, dt)
168
transientGeomCoeff=self._getTransientGeomCoeff(var),169

diffusionGeomCoeff=self._getDiffusionGeomCoeff(var),--> 170

buildExplicitIfOther=self._buildExplcitIfOther)171 172
self._buildCache(matrix, RHSvector)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/coupledBinaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
120

transientGeomCoeff=uncoupledTerm._getTransientGeomCoeff(tmpVar),
121

diffusionGeomCoeff=uncoupledTerm._getDiffusionGeomCoeff(tmpVar),-->
122
  buildExplicitIfOther=buildExplicitIfOther)123
 124 termMatrix += tmpMatrix
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/binaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
66
transientGeomCoeff=transientGeomCoeff, 67

diffusionGeomCoeff=diffusionGeomCoeff,---> 68

buildExplicitIfOther=buildExplicitIfOther) 69  70
matrix += tmpMatrix
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/unaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
97dt=dt,
98
transientGeomCoeff=transientGeomCoeff,---> 99
  diffusionGeomCoeff=diffusionGeomCoeff)
 100 elif buildExplicitIfOther:101 _, matrix,
RHSvector = self._buildMatrix(self.var,
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/abstractDiffusionTerm.pyc
in _