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 ]
>>>
