Dear Henri, For setting the values of CellVariables (e.g. initial conditions or variable coefficients), I recommend the following procedure, which avoids thinking about how the finite volume cells are organized in computer memory (and the associated headaches). Try this in a Jupyter Notebook or so.
I assume the follow imports. import numpy as np import fipy as fp Create your cylindrical (r,z) grid. Also, retrieve the 1D arrays which contain the r and z coordinates of each finite-volume cell. You may have noticed that all cells are stored in a 1D array, even for 2D and 3D grids. Note that the coordinates are the coordinates of the cell centers. gr = fp.CylindricalGrid2D(nr = 10 ,dr = 0.1, nz = 8, dz = 0.1) rr = gr.cellCenters[0] zz = gr.cellCenters[1] Create your relaxation CellVariable. And create a "Viewer" to see what you're doing. (The Viewer was repaired in a recent Fipy release, so please use the latest Fipy) tau1relaxs = fp.CellVariable(name = "Relaxation rate", mesh = gr ) vw = fp.Viewer(tau1relaxs) Fill all cells with the value 10.0. tau1relaxs[:] = 10. Now the magic begins. There is a syntax for conditionally indexing array elements. The following fills all cells situated at r < 0.5 with the value 2.0. Ask the Viewer object to generate a plot of the CellVariable. tau1relaxs[rr < 0.5] = 2. vw.plot() This can be elaborated further with element-wise 'and' (&). # make a disc radius 0.7 between z = 0.3 and z = 0.4 tau1relaxs[(rr < 0.7) & (zz > 0.3) & (zz < 0.4)] = 6. vw.plot() Hope this helps. Best regards, Martin On 10/11/2020 17:33, Henri Colaux wrote: > Dear Jon, dear Martinus, > > I really appreciate your replies! No worries about the delay. It is a > dificult period for all of us :( > > Thanks for letting me know that "zero-flux" (the condition I want) is > the implicit condition. I have no need to specify it specifically, I > just was not sure. > > I tried "CylindricalGrid2D" indeed, but I was not sure whether the > rotation axes was located at the left of the array, so I decided to > use "Grid2D" instead along with a "radius" variable so I am sure of > the position of the symmetry axis. Which brings me to my first > question: Here is one of the parameter arrays (here, the relaxation > time) I use in my MWE, before "ravel()" is apply: > > array([[ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.], > [ 2., 2., 2., 2., 2., 10., 10., 10.]]).ravel() > > Could you confirm that the symmetry axis is indeed located on the left > of this array ? > > In the MWE, I set domains 1 and 2 to have the same caracteristics for > simplicity. So, if I was correct of the way my system is programmed, I > should have for any time lines with identical values in the > polarisation variable "Pol" because there should be a cylindrical > symmetry. So it should give identical results to "CylindricalGrid1D". > Yet, if I sweep 50 times (for example), here is what I get : > > Pol_list[50] = > [[0.79 0.79 0.80 0.99 1.19 9.95 9.95 9.89] > [0.96 0.80 0.81 0.98 1.19 9.90 9.90 9.89] > [1.15 0.98 0.81 0.97 1.34 9.91 9.90 9.89] > [1.14 0.98 0.82 0.98 1.17 9.90 9.90 9.88] > [1.12 0.98 0.82 0.99 1.17 9.90 9.90 9.95] > [0.98 0.98 0.82 0.99 1.18 9.90 9.90 9.89] > [1.15 0.99 0.82 0.99 1.20 9.90 9.90 9.89] > [1.14 0.98 0.82 0.98 1.34 9.91 9.90 9.89] > [1.14 0.98 0.81 0.98 1.17 9.90 9.95 9.93] > [1.13 0.98 0.80 0.80 0.99 9.90 9.95 9.99]] > > Here, not only are my lines not identical, but the polarisation goes > above the equilibrium value of 1 inside the cylinder - the leftmost > values. This does not makes sens unless there is some transfer taking > place from the left of the system. Likewise, the rightmost values > should by the maximum, yet there are not. I tried with lower time > increments and the result stays the same. > > I have inclosed the updated MWE with no boundary condition and with > CylindricalGrid2D that gave the result above. Could you help my fix it > so it gives the same result as would an identical system with > CylindricalGrid1D ? > > Thanks in advance for your time if you have some :) > > PS: I am happy to hear there is now a « SphericalGrid1D » mesh which I > missed when starting my program (I have 3.3). > > Le mardi 10 novembre 2020 à 16:24:40 UTC+1, [email protected] a > écrit : > > Dear Henri, > > From your last message I gather that you would like to apply > zero-flux conditions on all boundaries, in a cylindrical (r,z) > geometry. > > In that case, I warmly recommend that you follow Jon's suggestion > and use the implicit no-flux boundary condition. Do not explicitly > specify a zero-gradient. We had some advection-diffusion test > cases here where explicitly specifying zero-flux conditions in > Fipy yielded less precise results (loss of mass conservation) than > the implicit no-flux. > > > https://github.com/mhvwerts/MasonWeaver-finite-diff/blob/master/finite-difference-solver-vs-Fipy-finite-volume-solver.pdf > > > Also, the CylindricalGrid2D is likely exactly what you need for > your problem geometry. We used it with success on a system for > which an analytic solution exists. > > > https://github.com/simulkade/FVTool/blob/master/Examples/External/Diffusion1DSpherical_Analytic-vs-FVTool-vs-Fipy/diffusion1Dspherical_analytic_vs_FVTool_vs_Fipy.pdf > > > Best wishes > Martin > > > P.S.1. Our example where we used CylindricalGrid2D was created > before SphericalGrid1D was implemented in Fipy > > P.S.2. For some it is a bit awkward to not specify boundary > conditions explicitly. It might be an idea to have a method > 'standard_zero_flux_boundaries()' or so, that simply does nothing, > but that one could put in one's code so that everything is > specified... > > > > On 10/11/2020 15:45, 'Guyer, Jonathan E. Dr. (Fed)' via fipy wrote: >> Henri - >> >> My apologies for the delay in answering. >> >> We have an existing mesh generator >> >> [`CylindricalGrid2D`](https://www.ctcms.nist.gov/fipy/fipy/generated/fipy.meshes.html#fipy.meshes.factoryMeshes.CylindricalGrid2D) >> that automatically takes care of calculating cylindrical >> symmetry. You only need to specify the r and z dimensions of the >> mesh and FiPy takes care of the wedge-shaped geometry; you do not >> need to transform your PDE into cylindrical form. >> >> Mathematically, what boundary conditions are you trying to achieve? >> >> I definitely wouldn’t use `NthOrderBoundaryCondition`. Those are >> only for higher-order PDEs like Cahn-Hilliard/Ginsburg-Landau. >> Even then, that’s an older way of doing things that we don’t >> recommend anymore. >> >> `Pol.faceGrad.constrain()` is the syntax we recommend, but you >> should be aware that as a finite volume code, FiPy has no-flux as >> the natural boundary condition. While no-flux is not always >> equivalent to zero normal gradient, for your diffusion equation, >> they are equivalent. There is no need to apply a constraint if >> no-flux is what you desire. >> >> - Jon >> >>> On Nov 3, 2020, at 12:49 PM, Henri Colaux <[email protected]> >>> wrote: >>> >>> >>> Dear devs, >>> >>> I am writing a code that simulates the diffusion of >>> magnetisation from proton to proton inside a solid. This process >>> behaves like a standard diffusion program (but with an >>> equilibrium term) so your code have been very helpful and I got >>> some great results so far, but I stambled into some problem >>> whose solution is above my knowledge. >>> >>> I am trying to simulate a « drill-core » sort of structure with >>> two different components - denoted « domains », one being a >>> small slice inside a much more abundant one, as pictured in the >>> following schematics : >>> >>> >>> I figured that this system is equivalent to producing a 2D array >>> using Grid2D as picture of the right of this figure, and >>> ensuring that there is no magnetisation transfer (i.e. no flux) >>> taking place at the border of the system. To account for the >>> rotation centre and the cylindrical symmetry, I have taken the >>> effect of the radius directly in the differential equation.Yet, >>> for some reasons, it seems like the boundery conditions is not >>> set correctly in my code, which gave unwanted behaviours. >>> >>> I have included a minimum working example of the code I am >>> trying to simulate. First and formost, could you let me know if >>> I do everything right for what I am trying to simulate ? >>> >>> Cheers ! >>> >>> Henri Colaux >>> >>> -- >>> To unsubscribe from this group, send email to >>> [email protected] >>> >>> View this message at https://list.nist.gov/fipy >>> --- >>> To unsubscribe from this group and stop receiving emails from >>> it, send an email to [email protected]. >>> <spindiff_v3_6_1_mwe.py> >> >> -- >> To unsubscribe from this group, send email to >> [email protected] >> >> View this message at https://list.nist.gov/fipy >> --- >> To unsubscribe from this group and stop receiving emails from it, >> send an email to [email protected]. > > -- To unsubscribe from this group, send email to [email protected] View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to [email protected].
