date:20161004

Re: how to set up data transfer between two adjacent nonuniform meshs

2016-10-04 Thread Zhekai Deng

Hi Ian,

Thank you very much for the example and pointing out the additional
reference. Those are very useful informations. This approach seems
intuitive and easy to incorporate into the code. I will try to incorporate
it into my own code.  Thanks again.

Best,

Zhekai

On Tue, Oct 4, 2016 at 5:23 AM, Campbell, Ian 
wrote:

> Hi Zhekai,
>
>
>
> There is.
>
>
>
> Try the following, where as you suggest, the interpolated value at the
> right-most face of mesh_1 is used. That (outward flowing) flux value is
> then copied to the inlet boundary condition (Neumann/flux) for mesh_2. A
> diffusion equation with a different coefficient for your solution variable
> phi is defined on the 2nd mesh – this should work if you change the
> equation further. Both meshes are non-uniform and of unit length.
>
>
>
> I don’t know from your question what type or magnitude boundary conditions
> you’re interested in on the other faces, but nonetheless, you can see in
> the figures produced from the code below that the outlet flux on the
> right-most face of mesh_1 is equal to the inlet flux on the left-most face
> of mesh_2. There’s likely a cleaner way to do this, by updating the value
> of the boundary condition for mesh_2 rather than re-defining it, but this
> is a start.
>
>
>
> An additional reference for you, here is Dr. Guyer’s suitably-named Gist:
> https://gist.github.com/guyer/bb199559c00f6047d466daa18554d83d
>
>
>
> Let us know if that works for you.
>
>
>
> Best,
>
>
>
> -  Ian
>
>
>
> phi_1_initial, phi_2_initial = 0.0, 0.0
>
> D1, D2 = 1.5, 1.75
>
> inlet_BC_value_init = 2.0
>
>
>
> dx_1 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
> 0.16221167, 0.10838638, 0.03806023]
>
> dx_2 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172,
> 0.16221167, 0.10838638, 0.03806023]
>
>
>
> mesh_1 = Grid1D(dx=dx_1, Lx=1.0)
>
> mesh_2 = Grid1D(dx=dx_2, Lx=1.0)
>
>
>
> phi_1 = CellVariable(mesh_1, value=phi_1_initial)
>
> phi_2 = CellVariable(mesh_2, value=phi_2_initial)
>
>
>
> eq_1 = TransientTerm() == ExplicitDiffusionTerm(coeff=D1)
>
> eq_2 = TransientTerm() == ExplicitDiffusionTerm(coeff=D2)
>
>
>
> # Boundary Conditions
>
> phi_1.faceGrad.constrain(0.5, where=mesh_1.facesLeft)
>
> phi_1.constrain(1.0, where=mesh_1.facesRight)
>
> phi_2.faceGrad.constrain(inlet_BC_value_init, where=mesh_2.facesLeft)
>
> phi_2.faceGrad.constrain(-0.5, where=mesh_2.facesRight)
>
>
>
> timeStepDuration = 0.9 * min(dx_1)**2 / (2 * max(D1, D2))
>
> steps = 100
>
> viewer = Viewer(vars=(phi_1, phi_2))
>
>
>
> for step in range(steps):
>
> eq_1.solve(var=phi_1, dt = timeStepDuration)
>
> inlet_BC_value = phi_1.faceValue[mesh_2.faceCenters == 1.0] #
> Acquire the data
>
> phi_2.faceGrad.constrain(inlet_BC_value, where=mesh_2.facesLeft)#
> Apply that data on the 2nd mesh
>
> eq_2.solve(var=phi_2, dt = timeStepDuration)
>
> viewer.plot()
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Zhekai Deng
> *Sent:* 04 October 2016 03:29
> *To:* fipy@nist.gov
> *Subject:* how to set up data transfer between two adjacent nonuniform
> meshs
>
>
>
> Hello All,
>
>
>
> I wonder is there any way to allow the data share on the two nonuniform
> mesh's interface ? To explain this, let's say I have two meshes, mesh_1 and
> mesh_2. Different equations govern the mesh_1 and mesh_2, and solution
> variable on both meshes is phi.  I would like the outlet (on the right
> face) on mesh_1 served as inlet(on the left face) on mesh_2.
>
>
>
> I tired to concatenate two mesh. If they are uniform, this works fine.
> However, if two does not share the exact the same face (for example, one is
> uniform, and another is nonuniform), there seems to be problem. .
>
>
>
> Thus, I wonder is there any way to "interpolate" the phi.facevalue on
> mesh_1 outlet face, and apply this interpolation value into the inlet of
> mesh_2?
>
>
>
>
>
> Best,
>
>
>
> Zhekai
>
> ___
> 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: sources errors in advection/diffusion problems, and one solution

2016-10-04 Thread James Pringle

Thank you.
I am busy for the week, but then will get a chance to think about this :-)

Cheers,
Jamie

On Tue, Oct 4, 2016 at 1:16 PM, Daniel Wheeler 
wrote:

> Hi James,
>
> Thanks for this example. I think I understand what is happening. Let's
> start with the code for solving the problem,
>
> 
> import fipy
> import numpy
>
> mesh = fipy.Grid1D(nx=10)
>
> var = fipy.CellVariable(mesh=mesh)
> var[:] = 2 * numpy.random.random(mesh.numberOfCells)
> var.constrain(1., where=mesh.facesRight)
> var.faceGrad.constrain(0, where=mesh.facesLeft)
>
> eqn = fipy.ConvectionTerm([-1]) == fipy.DiffusionTerm(1.0)
>
> for _ in range(100):
> eqn.sweep(var)
> print var
> ~~~
>
> Note that the "var.faceGrad.constrain" is required with convection
> terms to make sure that there is an outlet flow. This is a weird
> decision in many ways, but was made for backwards compatibility with
> previous FiPy versions, but this is not the main issue here.
>
> Notice that when the diffusion coefficient is small or large then this
> system solves in one sweep or close to one sweep (set it to 1000.0 or
> 0.0001 for example). The number of sweeps required to solve the system
> is at its worst when the diffiusion and convection coefficients are
> close in magnitude. So what is going on? Basically, the system is not
> fully implicit.
>
> The left hand boundary condition is dependent on the left hand cell
> value but gets added to the b vector in the case of non upwind
> convection terms. In the case of fully upwind convection terms this
> contribution is zero so as the diffusion coefficient goes to zero the
> CovectionTerm becomes fully upwind (or if a UpwindConvectionTerm is
> selected). As the diffusion coefficient gets large the contribution of
> the convection term decreases and thus the contribution from the
> explicit boundary condition is reduced.
>
> The discretization for the left most cell for a central difference
> scheme (we're using power law above, but it has some central
> difference contribution at D=1 and u=-1) is
>
>   D * phi0 - (u / 2 + D) * phi1 = -u * phiB / 2
>
> where phi0 is the value at the left most cell and phi1 is the cell to
> the right. In FiPy, phiB (the boundary value of phiB) is not
> calculated implicitly, but is explicitly set for the convection term
> boundary condition. Basically phiB is the previous value of phi0 in
> the FiPy code. That may not be the best way to deal with this. It may
> be that it needs to be explicit for stability reasons.
>
> A possible way to deal with this issue is to use terms for the outlet
> condition that are upwind at the outflow, but power law everywhere
> else.
>
> Anyway, I hope that helps with understanding what's going on even if
> the implementation is sub-optimal.
>
> Cheers,
>
> Daniel
>
> On Thu, 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

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

2016-10-04 Thread Daniel Wheeler

Hi James,

Thanks for this example. I think I understand what is happening. Let's
start with the code for solving the problem,

import fipy
import numpy

mesh = fipy.Grid1D(nx=10)

var = fipy.CellVariable(mesh=mesh)
var[:] = 2 * numpy.random.random(mesh.numberOfCells)
var.constrain(1., where=mesh.facesRight)
var.faceGrad.constrain(0, where=mesh.facesLeft)

eqn = fipy.ConvectionTerm([-1]) == fipy.DiffusionTerm(1.0)

for _ in range(100):
eqn.sweep(var)
print var
~~~

Note that the "var.faceGrad.constrain" is required with convection
terms to make sure that there is an outlet flow. This is a weird
decision in many ways, but was made for backwards compatibility with
previous FiPy versions, but this is not the main issue here.

Notice that when the diffusion coefficient is small or large then this
system solves in one sweep or close to one sweep (set it to 1000.0 or
0.0001 for example). The number of sweeps required to solve the system
is at its worst when the diffiusion and convection coefficients are
close in magnitude. So what is going on? Basically, the system is not
fully implicit.

The left hand boundary condition is dependent on the left hand cell
value but gets added to the b vector in the case of non upwind
convection terms. In the case of fully upwind convection terms this
contribution is zero so as the diffusion coefficient goes to zero the
CovectionTerm becomes fully upwind (or if a UpwindConvectionTerm is
selected). As the diffusion coefficient gets large the contribution of
the convection term decreases and thus the contribution from the
explicit boundary condition is reduced.

The discretization for the left most cell for a central difference
scheme (we're using power law above, but it has some central
difference contribution at D=1 and u=-1) is

  D * phi0 - (u / 2 + D) * phi1 = -u * phiB / 2

where phi0 is the value at the left most cell and phi1 is the cell to
the right. In FiPy, phiB (the boundary value of phiB) is not
calculated implicitly, but is explicitly set for the convection term
boundary condition. Basically phiB is the previous value of phi0 in
the FiPy code. That may not be the best way to deal with this. It may
be that it needs to be explicit for stability reasons.

A possible way to deal with this issue is to use terms for the outlet
condition that are upwind at the outflow, but power law everywhere
else.

Anyway, I hope that helps with understanding what's going on even if
the implementation is sub-optimal.

Cheers,

Daniel

On Thu, 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

Re: JOR solver to help in accelerating sweeps?

2016-10-04 Thread Daniel Wheeler

Hi Krishna,

Have you tried thoroughly profiling your code? You may want to check
how much of the time is being spent making copies of arrays and array
arithmetic as opposed to the linear solve before changing the solvers
for efficiency purposes. I'm not sure the linear solver matters for
reducing the number of sweeps. It may require a better non-linear
solver and using a single matrix. In the long run, it might be worth
contemplating some of the non-linear solvers in Trilinos that remove
the boundary between linear and non-linear solvers and also use matrix
free methods.

Sorry that I can't be more helpful.

On Thu, Sep 15, 2016 at 11:06 AM, Gopalakrishnan, Krishnakumar
 wrote:
>
> Hello,
>
>
>
> We have a specific question about how a solver helps in accelerating 
> convergence of a loosely coupled system of PDEs.
>
>
>
> Due to the fundamentally different properties/behaviour/geometry, we have a 
> loosely coupled system of 8 PDEs defined along 5 different meshes, with data 
> manually copied from certain regions of one mesh to suitable regions of 
> another after sweeping each PDE.
>
>
>
> Quite unsurprisingly, the system of equations needs a lot of sweeps (80 to 
> 100) within each time-step to converge, and consequently the simulation model 
> isn’t useable. Although, the trilinos based parallel solvers are correctly 
> set-up and working correctly for the example problems shipped with FiPy, the 
> mpirun command complains when invoked on our (complicated script). We assume 
> that it is due to the complicated setup of meshes and the classes/methods and 
> objects that we have instantiated for our application. The serial code runs 
> just fine.
>
>
>
> So, in order to speed up our simulation, we looked to implement dynamic 
> relaxation factors for the sweeps. This has only a limited success so far.
>
>
>
> Next, I wonder whether including a Successive over-relaxation solver will 
> help in this case ? There are 8 PDEs to be swept continuously until the 
> residue of all of them gets lowered below a tolerance.  Does the Fipy’s JOR 
> solver (wrapper for Pysparse Jacobi/over-relaxation solvers)  help in this 
> case ?  is there any other way to speed up our system ?
>
>
>
>
>
> Best Regards
>
>
>
> Krishna
>
>
>
>
>
>
>
>
>
>
> ___
> 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: how to set up data transfer between two adjacent nonuniform meshs

2016-10-04 Thread Campbell, Ian

Hi Zhekai,

There is.

Try the following, where as you suggest, the interpolated value at the 
right-most face of mesh_1 is used. That (outward flowing) flux value is then 
copied to the inlet boundary condition (Neumann/flux) for mesh_2. A diffusion 
equation with a different coefficient for your solution variable phi is defined 
on the 2nd mesh – this should work if you change the equation further. Both 
meshes are non-uniform and of unit length.

I don’t know from your question what type or magnitude boundary conditions 
you’re interested in on the other faces, but nonetheless, you can see in the 
figures produced from the code below that the outlet flux on the right-most 
face of mesh_1 is equal to the inlet flux on the left-most face of mesh_2. 
There’s likely a cleaner way to do this, by updating the value of the boundary 
condition for mesh_2 rather than re-defining it, but this is a start.

An additional reference for you, here is Dr. Guyer’s suitably-named Gist: 
https://gist.github.com/guyer/bb199559c00f6047d466daa18554d83d

Let us know if that works for you.

Best,


-  Ian

phi_1_initial, phi_2_initial = 0.0, 0.0
D1, D2 = 1.5, 1.75
inlet_BC_value_init = 2.0

dx_1 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172, 0.16221167, 
0.10838638, 0.03806023]
dx_2 = [0.03806023, 0.10838638, 0.16221167, 0.19134172, 0.19134172, 0.16221167, 
0.10838638, 0.03806023]

mesh_1 = Grid1D(dx=dx_1, Lx=1.0)
mesh_2 = Grid1D(dx=dx_2, Lx=1.0)

phi_1 = CellVariable(mesh_1, value=phi_1_initial)
phi_2 = CellVariable(mesh_2, value=phi_2_initial)

eq_1 = TransientTerm() == ExplicitDiffusionTerm(coeff=D1)
eq_2 = TransientTerm() == ExplicitDiffusionTerm(coeff=D2)

# Boundary Conditions
phi_1.faceGrad.constrain(0.5, where=mesh_1.facesLeft)
phi_1.constrain(1.0, where=mesh_1.facesRight)
phi_2.faceGrad.constrain(inlet_BC_value_init, where=mesh_2.facesLeft)
phi_2.faceGrad.constrain(-0.5, where=mesh_2.facesRight)

timeStepDuration = 0.9 * min(dx_1)**2 / (2 * max(D1, D2))
steps = 100
viewer = Viewer(vars=(phi_1, phi_2))

for step in range(steps):
eq_1.solve(var=phi_1, dt = timeStepDuration)
inlet_BC_value = phi_1.faceValue[mesh_2.faceCenters == 1.0] # 
Acquire the data
phi_2.faceGrad.constrain(inlet_BC_value, where=mesh_2.facesLeft)# Apply 
that data on the 2nd mesh
eq_2.solve(var=phi_2, dt = timeStepDuration)
viewer.plot()

From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Zhekai 
Deng
Sent: 04 October 2016 03:29
To: fipy@nist.gov
Subject: how to set up data transfer between two adjacent nonuniform meshs

Hello All,

I wonder is there any way to allow the data share on the two nonuniform mesh's 
interface ? To explain this, let's say I have two meshes, mesh_1 and mesh_2. 
Different equations govern the mesh_1 and mesh_2, and solution variable on both 
meshes is phi.  I would like the outlet (on the right face) on mesh_1 served as 
inlet(on the left face) on mesh_2.

I tired to concatenate two mesh. If they are uniform, this works fine. However, 
if two does not share the exact the same face (for example, one is uniform, and 
another is nonuniform), there seems to be problem. .

Thus, I wonder is there any way to "interpolate" the phi.facevalue on mesh_1 
outlet face, and apply this interpolation value into the inlet of mesh_2?


Best,

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

Re: how to set up data transfer between two adjacent nonuniform meshs

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

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

Re: JOR solver to help in accelerating sweeps?

RE: how to set up data transfer between two adjacent nonuniform meshs

5 matches

Site Navigation

Mail list logo

Footer information