Martin - That's interesting that `u.value.reshape(Nz,Ny,Nx)` gives you a view on the CellVariable's value. I'd say it's not 100% reliable that it always will, but it probably will. The result is certainly easier to read than juggling around with `ravel_multi_index`.
- Jon > On Apr 8, 2019, at 12:01 PM, Martinus WERTS <[email protected]> > wrote: > > Dear Dario, > > I am currently using Fipy for a slightly different diffusion-reaction > problem. > > I use a CellVariable on the same grid for the (spatially-dependent) rate > constant. One can set its values in a "numpy-vectorized" fashion by accessing > the "value" property of the CellVariable. > > By using a Boolean expression for indexing the "value" vector one can address > individual elements. > > Here are bits of code that illustrate the approach. It sets the "capture > rate" to Kcap inside a sphere of radius R. Outside the sphere, "capture rate" > is zero. > > # create grid (one octant, using symmetry) > mesh = UniformGrid3D(nx=Nx, ny=Ny, nz=Nz, > dx=sx/Nx, dy=sy/Ny, dz=sz/Nz, > origin=[[0.],[0.],[0.]] ) > X,Y,Z = mesh.cellCenters > > # variables > u = CellVariable(name="particles", > mesh=mesh, > value=0.) > k = CellVariable(name="capture rate", > mesh=mesh, > value=0.) > > # equation > eq = TransientTerm(var = u) == DiffusionTerm(coeff = Dcoeff, var = u) - > k*u > > # I skip the boundary conditions here to keep this message compact > > # initial conditions for u; setting value for rate constant > u.value=Ceq > k.value[(X**2 + Y**2 + Z**2< Rcap**2)] = Kcap > > > I also found that it is possible to access the elements on this regular 3D > grid via a "reshaped" version of the array > > uvol = u.value.reshape(Nz,Ny,Nx) > > (...) > > uvol[:,0,:] += np.sum(vvol, axis=1) > > This approach, however, requires that one does the book-keeping of the > coordinate system oneself, and relies on the specific way that the data is > organized by FiPy. > > Best wishes, > Martin > > > > > > On 08/04/2019 17:30, Dario Panada wrote: >> Hello, >> >> Many thanks for your reply. >> >> Yes, let me provide a bit more context. >> >> I two initial numpy grids (n*n*n) where each value corresponds to a >> source/sink. Eg: Given my source grid and coordinates (1,2,3) having value >> 5, I want to set such value as a source in FiPy. Currently I am dong that >> by, for each such coordinate, finding the corresponding ravelled >> index and setting it in _array, as snippet in my previous message; >> >> i = np.ravel_multi_index([coordinate[0], coordinate[1], coordinate[2]], >> (20, 20, 20)) >> >> sourceGrid._array[i] = sourceRate >> sinkGrid._array[i] = sinkRate >> >> I suppose I could build the entire vector of sources before and then doing a >> single assignment to _array, but again you correctly mentioned >> that relying on that is bad practice. >> >> You mention: >> >> sourceGrid[..., i] = sourceRate >> >> Can I just please confirm what data type sourceGrid is? In the context of >> defining the equation >> >> eq = TransientTerm() == DiffusionTerm(coeff=D) + sourceGrid - sinkGrid >> >> Can sourceGrid/sinkGrid just be numpy arrays or even simple python lists? I >> was under the impression they had to be CellVariable objects but could be >> wrong. >> >> Kind Regards, >> Dario >> >> >> >> On Mon, Apr 8, 2019 at 3:04 PM Guyer, Jonathan E. Dr. (Fed) via fipy >> <[email protected]> wrote: >> Iterating over a mesh with a Python `for` loop is, as you've found, an >> incredibly inefficient way to do things. FiPy, like numpy it relies on, is >> intended to be used with vectorized operations. >> >> As far as your approach, things that start with `_` in Python are internal >> implementation details and you should not depend on them. If you find cases >> in FiPy that absolutely require that you access `some._propertyName` or >> `some._methodName()`, then please file an issue explaining your need so that >> we can provide a public interface. >> >> In this case, there's no need to access `_array`. Just write >> >> >>> sourceGrid[..., i] = sourceRate >> >>> sinkGrid[..., i] = sinkRate >> >> HOWEVER, you are still relying on an internal implementation detail, >> specifically that a Grid is an array with known fastest and slowest varying >> axes. We have in FiPy's history switched from Fortran to C ordering, and we >> might conceivably switch back at some point. Further, if you every wanted to >> run your script on an unstructured mesh, then your system wouldn't work at >> all. >> >> `sourceCoords` presumably comes from some definition in terms of geometry >> rather than discrete mesh indices. FiPy is intended to work best with those >> geometric descriptions. If you described where `sourceCoords` came from, we >> could help with a more FiPyish way to get what you want. >> >> >> >> > On Apr 6, 2019, at 9:02 AM, Dario Panada <[email protected]> wrote: >> > >> > Good Afternoon (Morning) to all, >> > >> > I have an equation of type >> > >> > eq = TransientTerm() == DiffusionTerm(coeff=D) + sourceGrid - sinkGrid >> > >> > Where sourceGrid and sinkGrid are derived from values in a 3D grid. >> > >> > Is there any downside to declaring the grids as >> > >> > sourceGrid = CellVariable(name="source", mesh=Grid3D(dx=1, dy=1, dz=1, >> > nx=20, ny=20, nz=20)) >> > sinkGrid = CellVariable(name="sink", mesh=Grid3D(dx=1, dy=1, dz=1, >> > nx=20, ny=20, nz=20)) >> > >> > And populating them as: >> > >> > i = np.ravel_multi_index([coordinate[0], coordinate[1], coordinate[2]], >> > (20, 20, 20)) >> > >> > sourceGrid._array[i] = sourceRate >> > sinkGrid._array[i] = sinkRate >> > >> > The original procedure I was using (given below) where I called .setValue >> > on each coordinate is extremely time-consuming. Results seem to align, but >> > of course this doesn't mean it's right... Or if it's wrong, any specific >> > advice on how to achieve this? >> > >> > Thanks, >> > Dario >> > >> > Original approach (and there is an equivalent function for setupSinkGrid): >> > def setupSourceGrid_(self, sourceCoords, mesh): >> > sourceGrid = CellVariable(name="source", mesh=mesh, value=0) >> > sourceGrid.setValue(0.) >> > for pos, v in sourceCoords.iteritems(): >> > tmpGrid = [False for _ in range(8000)] >> > i = np.ravel_multi_index([pos[0], pos[1], pos[2]], (20, 20, >> > 20)) >> > tmpGrid[i] = True >> > sourceGrid.setValue(v, where=tmpGrid) >> > >> > _______________________________________________ >> > fipy mailing list >> > [email protected] >> > http://www.ctcms.nist.gov/fipy >> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >> >> >> _______________________________________________ >> fipy mailing list >> [email protected] >> http://www.ctcms.nist.gov/fipy >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >> >> >> _______________________________________________ >> fipy mailing list >> >> [email protected] >> http://www.ctcms.nist.gov/fipy >> >> [ NIST internal ONLY: >> https://email.nist.gov/mailman/listinfo/fipy >> ] >> > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
