RE: Applying Bounding Limits to PDE Solution Variables
Hi Daniel, Thank you very much for your prompt reply. I'm manually trying to implement the bound-checking feature immediately following a sweep, as you've suggested. As an initial step towards this, I need to manually compute my own residual and compare to that computed via '.sweep' in FiPy. I've tried the mean-square and RMS values, but they're both orders of magnitude away from the residual returned by '.sweep'. As such, can you please tell me the formula used by '.sweep' to compute the residual, so that I may use the same one? With best regards, - Ian -Original Message- From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 02 August 2016 15:57 To: Multiple recipients of list Subject: Re: Applying Bounding Limits to PDE Solution Variables Hi Ian, Sorry, there isn't any way to do this as an implicit command in FiPy, whether as a simple reset during sweeping or as a fully implicit constraint in the equations. Of course, you can always try these techniques explicitly in the Python code or in the design of the equations. I agree that It would be nice to have a one sided constraint syntactically in FiPy, but I'm not sure how best to make that work internally. Cheers, Daniel On Tue, Aug 2, 2016 at 10:46 AM, Campbell, Ian wrote: > Hi All, > > > > I would like to prevent the CellVariables I am solving for from ever > reaching physically-meaningless and incorrect values, such as > concentrations becoming negative. Is there a way to do this? > > > > The most natural place I could think of it being is an argument in > CellVariable declaration, such as “phi = CellVariable(mesh = mesh2d, > value = phi_0, upperlimit = 10**7, lowerlimit = 0)”. However, I > haven’t found such an option there or in any of the FiPy documentation. > > > > With best regards, > > > > - Ian > > > ___ > 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 ] ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Applying Bounding Limits to PDE Solution Variables
Hi Ian, Sorry, there isn't any way to do this as an implicit command in FiPy, whether as a simple reset during sweeping or as a fully implicit constraint in the equations. Of course, you can always try these techniques explicitly in the Python code or in the design of the equations. I agree that It would be nice to have a one sided constraint syntactically in FiPy, but I'm not sure how best to make that work internally. Cheers, Daniel On Tue, Aug 2, 2016 at 10:46 AM, Campbell, Ian wrote: > Hi All, > > > > I would like to prevent the CellVariables I am solving for from ever > reaching physically-meaningless and incorrect values, such as concentrations > becoming negative. Is there a way to do this? > > > > The most natural place I could think of it being is an argument in > CellVariable declaration, such as “phi = CellVariable(mesh = mesh2d, value = > phi_0, upperlimit = 10**7, lowerlimit = 0)”. However, I haven’t found such > an option there or in any of the FiPy documentation. > > > > With best regards, > > > > - Ian > > > ___ > 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 ]
Applying Bounding Limits to PDE Solution Variables
Hi All, I would like to prevent the CellVariables I am solving for from ever reaching physically-meaningless and incorrect values, such as concentrations becoming negative. Is there a way to do this? The most natural place I could think of it being is an argument in CellVariable declaration, such as "phi = CellVariable(mesh = mesh2d, value = phi_0, upperlimit = 10**7, lowerlimit = 0)". However, I haven't found such an option there or in any of the FiPy documentation. With best regards, - Ian ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: reaction-advection-diffusion; problems making the solution converge
Hi Nils, Could you possible post your code so I can play around with it? I can't think offhand what could be the issue or exactly what the error is that you're calculating, it seems like the normalized residual or something. This is a linear equation so I'm not entirely sure why the non-linear residual wouldn't be very close to zero if the linear solver is converging. Probably best just to see the code rather than speculate. Cheers, Daniel On Mon, Aug 1, 2016 at 7:50 AM, Nils Becker wrote: > hi, > > i scanned the list and the documentation but have not been able to make > this work properly yet: > > i have a 2D pde of the form > > \[ > > \partial_t \rho = \nabla \cdot (\rho \nabla \phi) + D \nabla^2 \rho + > \lambda \rho > > \] > > where phi is a potential field and lambda is a source field, both are > constant in time and smoothly varying in space. > > i have been solving this with fipy, using > > DiffusionTermCorrection(coeff=[D,]) for the diffusion > PowerLawConvectionTerm(coeff=v ) where v = -phi.grad for the advection > ImplicitSourceTerm(coeff=lambda) for the source > > the boundary conditions are either no-flux or dirichlet. i think > advection is dominant, although i can't say i have actually tried to > quantify that. > > i am testing the solution by subtracting the rhs and lhs of the equation > and comparing that to \partial_t \rho in magnitude. i am having trouble > making this error smaller than around 10%: at points where the density > grows strongly, the solution deviates. > > at the moment i am using a loop like this: > > for (st, stnext) in zip(Solve.sample_t, Solve.sample_t[1:]): > for tt in nx.arange(st, stnext, Solve.step): > rho.updateOld() > res = 1. > while res > Solve.tolerance: > res = eq.sweep(var=rho, dt=Solve.step) > > i've played around with making the mesh smaller -- did not seem to make > any difference; reducing the tolerance -- some effect, and reducing the > temporal step size, which did not help much. > > should i be doing something differently? > > thanks for any help! > > n. > ___ > 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: error in gradient even with orthogonal grid
Thank you -- since my problem solves for a stream function, this was a very
helpful tip, and has made my solutions appear much better.
Jamie
On Mon, Aug 1, 2016 at 4:49 PM, Daniel Wheeler
wrote:
> Hi Jamie,
>
> Sorry for the delayed reply.
>
> On Thu, Jul 21, 2016 at 11:31 AM, James Pringle wrote:
> >
> > However, I am still getting large errors in estimates of the gradient
> even
> > when the mesh is nearly perfectly orthogonal. These errors DO NOT become
> > smaller as resolution is increased. I have two question:
> >
> > Are these errors to be expected in unstructured meshes?
>
> Cell centered FV has two issues with unstructured meshes. You know
> about the orthogonality issue, but there is a second issue. The
> docstring from the following code might help explain,
>
>
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_usnistgov_fipy_blob_develop_fipy_variables_cellVariable.py-23L257&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=RM15FSI_eegFR5VfTb9il1PZ2B1_yaQ6eKmZHq-8Hx0&e=
>
> See, `\frac{1}{V_P} \sum_f \vec{n} \phi_f A_f`. This is the
> discretized gradient calculation. In this formulation using the
> triangular mesh, the calculated \phi_f is not the correct average for
> the face as the line between the cell centers doesn't cross the face
> in the middle of the face. I think this is called non-conjunctionality
> in finite volume speak. See equation 9.8 in the following link for
> further explanation (we should probably implement this approach).
>
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__bit.ly_2aqRiHt&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=EJ34X0r43RF50EgJ2p_clCS-H5yY8nSTitKGqtbINts&e=
>
> I believe that this is the source of the error in the gradient
> calculation. It's a consistent error independent of the mesh
> refinement I think (or maybe first order). An alternative gradient
> calculation is the least squares gradient, see
>
>
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_usnistgov_fipy_blob_develop_fipy_variables_cellVariable.py-23L270&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=l6F4maMJnB3bfM6--FXXLzrUKyaH1i7KHuo1atLYJqM&e=
>
> Using that
>
>
> import fipy as fp
> nx = 25
> dx = 1. / nx
> mesh = fp.Tri2D(nx=nx, ny=nx, dx=dx, dy=dx)
> psi = fp.CellVariable(mesh=mesh)
> psi.constrain(1.0, where=mesh.facesLeft)
> psi.faceGrad.constrain(1.0, where=mesh.facesRight)
> eq = fp.DiffusionTerm().solve(psi)
> print psi.leastSquaresGrad[:, :100]
>
>
> seems to give lots of 1s and 0s which seems correct. The cells next to
> the boundary are not correct as they aren't using the boundary
> conditions in the calculation. I'm not sure about the order of
> accuracy of the least squares or why it happens to work for this mesh.
>
> > Is it acceptable to average the gradient over space to reduce the error?
>
> I don't think so unless it can be shown to improve the accuracy with
> analysis.
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=A7D0acjE011uAk-Vl-5OhrTNVdH9kgmq0rshSN6tbsg&e=
> [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=V8y_Fr-VcNq0fTtFnfiyabj8MXL8mmGa5LFCMv0Sx_Q&e=
> ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
