[Dune] Periodicity for unstructured grids

[Matteo Semplice]

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