[Dune] Periodicity for unstructured grids

Matteo Semplice matteo.semplice at unito.it
Tue Dec 13 08:01:33 CET 2016 

Dear Martin,
     thanks a lot for the clarification (and for unveiling the 
dist=geoIn.global( xIn ) - geoOut.global( xOut )  trick!)

Cheers,
     Matteo

On 12/12/2016 20:31, Martin Nolte wrote:
> Hi Matteo,
>
> maybe I did not make myself clear: There are two suggestions on the 
> implementation of periodicity (see my first answer on this thread). 
> The documentation always refers to the only officially supported 
> approach of discretizing the n-torus (Approach), which is implemented 
> by YaspGrid only. AlbertaGrid, ALUGrid and SPGrid implement Approach 
> 2, which I usually refer to as "support for periodic boundary 
> conditions" to emphasize the difference to modeling the actual 
> manifold. The latter approach was never officially discussed and can, 
> hence, not be part of the documentation.
>
> It might be confusing that the majority of grids supporting periodic 
> boundary conditions using an unofficial approach. However, as you also 
> noticed, there has been little interest in the topic of periodicity 
> within the DUNE community, so the developer meeting never got around 
> to a proper discussion.
>
> To make this clear: There is only one officially supported approach to 
> periodicity: The one implemented in YaspGrid. The other three grids I 
> mentioned above only support an unofficial interface extension that 
> should be considered exprimental, though it is quite stable for the 
> last 7 years.
>
> Best,
>
> Martin
>
> On 12/12/2016 07:22 PM, Matteo Semplice wrote:
>> Hi,
>>
>>     maybe it's me being dumb, but
>>
>> https://www.dune-project.org/doxygen/2.4.1/classDune_1_1Intersection.html 
>>
>>
>> still says that on an intersection located on a periodic boundary both
>> neighbor() and boundary() return true and outside() returns a 
>> "Ghost-/Overlap
>> cell (with transformed geometry)".
>>
>> I understand the previous statement as follows. Say that I have a 1d 
>> periodic
>> grid for the interval [0,1], with elements of size h_i for 
>> i=0,...N-1. When I
>> ask for outside() from the interface located at x=0, I understant 
>> that I will
>> get a cell of size h_{N-1} but with center located at x<0. However, 
>> if I apply a
>> mapper to this outside() element, the mapper tells me that the dofs 
>> of the
>> outside element are located at the same place where the dofs of the 
>> real cell
>> number N-1 are stored. In this sense I was saying earlier that user 
>> code (for
>> finite volume explicit methods) does not even realize that it is 
>> operating on a
>> boundary: it can compute polynomial reconstructions, fluxes, etc as 
>> it would do
>> on an inner intersection. When I last tested, this part on the mapper 
>> did work
>> as I was expecting, but the previous one about the geometry did not.
>>
>> On the other hand Martin is now telling that in order to distinguish 
>> a periodic
>> boudary I shouldnot follow the table in the abovementioned docs, but 
>> instead
>> compute "dist" as suggested below and compare it with zero, and (I 
>> guess) make
>> the geometry transformation by myself if needed.
>>
>> I did not check the docs for the proposed 2.5 release, but may I 
>> suggest to
>> clarify this point (unless it has already been done, of course)?
>>
>> Best regards,
>>     Matteo
>>
>> On 12/12/16 16:01, Aleksejs Fomins wrote:
>>> Hey Martin,
>>>
>>> Thank you very much, that's exactly what I was looking for :)
>>>
>>>
>>> On 12.12.2016 15:00, Martin Nolte wrote:
>>>> Hi Aleksejs,
>>>>
>>>> for the "periodic boundary condition" approach, the periodic 
>>>> boundaries are
>>>> marked as boundaries (e.g, AlbertaGrid, ALUGrid, SPGrid), but 
>>>> report neighbor
>>>> to be true. They will return the "glued" element on the other side 
>>>> of the
>>>> domain as their outside entity.
>>>>
>>>> You can, then, obtain the distance between the two boundaries by 
>>>> the following
>>>> expression:
>>>>
>>>> auto dist = geoIn.global( xIn ) - geoOut.global( xOut )
>>>>
>>>> where
>>>> - geoIn = intersection.inside().geometry()
>>>> - geoOut = intersection.outside().geometry()
>>>> - xIn = intersection.geometryInInside().global( x )
>>>> - xOut = intersection.geometryInOutside().global( x )
>>>>
>>>> For non-periodic intersection, dist is zero.
>>>>
>>>> Notice: There is one little caveat: We usually require boundary 
>>>> intersections
>>>> to be comforming. For periodic boundaries, I always ignored this 
>>>> requirement.
>>>>
>>>> Best,
>>>>
>>>> Martin
>>>>
>>>> On 12/12/2016 10:10 AM, Aleksejs Fomins wrote:
>>>>> Dear Matteo,
>>>>>
>>>>> Thank you for your answer. While I see the intuitive motivation of 
>>>>> having
>>>>> translated periodic ghost elements, I would currently prefer to 
>>>>> keep the
>>>>> periodic neighbor there where it is and not translate it. In 
>>>>> Floquet's
>>>>> (Bloch's) theory, applicable to i.e. electromagnetics, the field 
>>>>> must obtain
>>>>> a phase factor depending on the distance between periodic 
>>>>> boundaries, and
>>>>> thus the user code must be aware of periodic neighbors and treat them
>>>>> differently from interior and ghost neighbors. In effect, creating
>>>>> translated elements increases grid complexity and storage, but I 
>>>>> do not see
>>>>> yet how it would simplify the user code. When calculating the 
>>>>> coupling
>>>>> matrices, integrated over the periodic boundary, the contributions 
>>>>> from each
>>>>> neighbor can be evaluated using the same local coordinate, given 
>>>>> that the
>>>>> elements are adjusted by the grid to exactly match the orientation.
>>>>>
>>>>> I have a question: How does one denote a periodic boundary in 
>>>>> Dune? As far
>>>>> as I am aware, the partition types are
>>>>>
>>>>> border,
>>>>> interior,
>>>>> front,
>>>>> overlap and
>>>>> ghost
>>>>>
>>>>> Is there some other variable within the grid that is responsible 
>>>>> for the
>>>>> periodicity?
>>>>>
>>>>> Best regards,
>>>>> Aleksejs
>>>>>
>>>>>
>>>>> On 12.12.2016 08:42, Matteo Semplice wrote:
>>>>>> Dear everybody,
>>>>>>
>>>>>>      a few years ago I had raised the issue that the behaviour of 
>>>>>> dune was
>>>>>> in fact different from what announced in the documentation. In 
>>>>>> particular I
>>>>>> read the documentation as saying "if you ask for the outer 
>>>>>> element from an
>>>>>> interface on a periodic boundary, you get a fake element with the
>>>>>> periodically transformed geometry but the id of the real 
>>>>>> element". This
>>>>>> would allow to make no change at all in the user code for periodic
>>>>>> boundaries, since the user would not see at all the boundary. 
>>>>>> Unfortunately
>>>>>> most grid implementations at the time did not comply to this 
>>>>>> statement and
>>>>>> the discussion on the list ended with "we will change the 
>>>>>> documentation",
>>>>>> which was entirely fair since I was not volunteering to do the 
>>>>>> job! So, on
>>>>>> the one hand, this is just to say that it's not that I lost 
>>>>>> interest, I
>>>>>> just stopped asking (since again I do not feel I have the ability 
>>>>>> to do the
>>>>>> job).
>>>>>>
>>>>>> On the other hand, at the time I had written a small code that 
>>>>>> can be
>>>>>> compiled with the usual "GRIDTYPE=XX GRIDDIM=NN make" trick and 
>>>>>> that checks
>>>>>> the behaviour of dune on periodic interfaces for that grid 
>>>>>> implementation.
>>>>>> Needless to say that I would be happy to share it, if it can be 
>>>>>> of interest
>>>>>> to anyone working on periodic boundaries.
>>>>>>
>>>>>> Best regards,
>>>>>>      Matteo
>>>>>>
>>>>>> On 09/12/16 12:53, Martin Nolte wrote:
>>>>>>> Hi Aleksejs,
>>>>>>>
>>>>>>> there has been quite some discussion on periodic boundary 
>>>>>>> conditions.
>>>>>>> Currently, there are two approaches:
>>>>>>>
>>>>>>> Approach 1:
>>>>>>>
>>>>>>> Consider the flat torus T^n = R^n / Z^n. In this case, you can 
>>>>>>> condider your
>>>>>>> computational domain to be R^n itself and theoretically imagine 
>>>>>>> a grid with
>>>>>>> infinitely many elements such that for each (geometrical) 
>>>>>>> element E and
>>>>>>> each i
>>>>>>> in Z^n, E + i is also an element of the grid. Factoring out Z^n, 
>>>>>>> you end up
>>>>>>> with a grid for T^n.
>>>>>>>
>>>>>>> Now, imagine you also have infinitely many processors and assign 
>>>>>>> to each cube
>>>>>>> [0,1]^n to one processor. This case can (theoretically) be 
>>>>>>> handled by the
>>>>>>> parallel interface. As each process only has a shifted copy of 
>>>>>>> the grid, we
>>>>>>> periodicity allows us to consider only one representative of the 
>>>>>>> parallel
>>>>>>> processes, i.e., 1 process (which then communicates with 
>>>>>>> itself). I hope
>>>>>>> it is
>>>>>>> clear how this can be extended to multiple "real" processors.
>>>>>>>
>>>>>>> This approach is implemented in YaspGrid.
>>>>>>>
>>>>>>> There is one little tweak to notice: As far as I know, the 
>>>>>>> identified
>>>>>>> elements
>>>>>>> are considered copies of each other. This might lead to problems 
>>>>>>> with the
>>>>>>> coordinate function, which is not periodic at all. I'm not sure 
>>>>>>> whether this
>>>>>>> question was thoroughly discussed, though.
>>>>>>>
>>>>>>>
>>>>>>> Approach 2:
>>>>>>>
>>>>>>> Let us consider the flat torus to be constructed by glueing 
>>>>>>> together opposite
>>>>>>> sides of the unit cube. In other words, we consider elements on 
>>>>>>> the right
>>>>>>> boundary to have an intersection with elements on the left 
>>>>>>> boundary (and vice
>>>>>>> versa). But that's it, the subentities are not identified.
>>>>>>>
>>>>>>> In this case, parallelization is also clear. For each 
>>>>>>> intersection, the
>>>>>>> neighboring element should exist as a ghost.
>>>>>>>
>>>>>>>
>>>>>>> A few years ago, there was quite some discussion on this, but 
>>>>>>> people seem to
>>>>>>> have lost interest. I'm not sure whether there is a "right and 
>>>>>>> wrong".
>>>>>>> Personally, I consider both approaches viable, depending on your 
>>>>>>> goal.
>>>>>>>
>>>>>>> Best,
>>>>>>>
>>>>>>> Martin
>>>>>>>
>>>>>>> PS: You can find some more information in the appendix of my phd 
>>>>>>> thesis.
>>>>>>>
>>>>>>>
>>>>>>> On 12/09/2016 11:52 AM, Aleksejs Fomins wrote:
>>>>>>>> Dear Dune,
>>>>>>>>
>>>>>>>> I need to implement periodic boundary conditions into
>>>>>>>> dune-curvilineargrid now.
>>>>>>>>
>>>>>>>> Would you be so kind to tell me if the current dune-grid interface
>>>>>>>> envisions periodic boundaries and associated ghost elements, 
>>>>>>>> and how.
>>>>>>>>
>>>>>>>> The questions of interest are: partiton type, behavior of 
>>>>>>>> intersection,
>>>>>>>> behavior of DataHandle
>>>>>>>>
>>>>>>>> Best regards,
>>>>>>>> Aleksejs
>>>>>>>>
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> Dune mailing list
>>>>>>>> Dune at dune-project.org
>>>>>>>> http://lists.dune-project.org/mailman/listinfo/dune
>>>>>>>>
>>>>>> _______________________________________________
>>>>>> Dune mailing list
>>>>>> Dune at dune-project.org
>>>>>> http://lists.dune-project.org/mailman/listinfo/dune
>>>>>
>>>>> _______________________________________________
>>>>> Dune mailing list
>>>>> Dune at dune-project.org
>>>>> http://lists.dune-project.org/mailman/listinfo/dune
>>>>>
>>>
>>> _______________________________________________
>>> Dune mailing list
>>> Dune at dune-project.org
>>> http://lists.dune-project.org/mailman/listinfo/dune
>>
>>
>> _______________________________________________
>> Dune mailing list
>> Dune at dune-project.org
>> http://lists.dune-project.org/mailman/listinfo/dune
>

More information about the Dune mailing list