Doubly-periodic BC syntax

Hi,

I may be wrong, but I don't think it will work out of the box with two periodic directions in one PeriodicBC object? I know there is nothing in there right now that handles those corners. I used to work with two periodic directions using two PeriodicBC objects, but you cannot do that anymore either.

The handling of periodic boundaries are about to change. If you are in a hurry you can try this branch here: lp:/~mikael-mortensen/dolfin/periodic_to_dofmap, where the periodic boundaries are incorporated directly into the dofmaps. It works with two and three periodic directions as well, but the code is under review and not mature nor supported yet.

Best regards

Mikael

Den Nov 15, 2012 kl. 3:11 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214267
>
> Description changed to:
> I am trying to solve a problem on a doubly-periodic domain. However,
> the answer is clearly wrong, as it has differing values at "equivalent"
> nodes at the corners. This strongly suggests that the 4 corners are not
> being contracted to a single degree of freedom. I have no idea whether
> my program is wrong, or whether this is a bug, so I'm asking whether my
> code below is correct.
>
> from dolfin import *
> set_log_level(ERROR)
> mesh = UnitSquare(1, 1)
> E = FunctionSpace(mesh, 'CG', 2)
> S = FunctionSpace(mesh, 'BDM', 1)
>
> class PeriodicBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and ( (x[0]<DOLFIN_EPS and x[1]<(1-DOLFIN_EPS)) or (x[1] < DOLFIN_EPS and x[0]<(1-DOLFIN_EPS)))
> def map(self, x, y):
> if (x[0] > (1 - DOLFIN_EPS)) and (x[1] > (1 - DOLFIN_EPS)):
> y[0] = x[0] - 1
> y[1] = x[1] - 1
> elif x[0] > (1 - DOLFIN_EPS):
> y[0] = x[0] - 1
> y[1] = x[1]
> else:
> y[0] = x[0]
> y[1] = x[1] - 1
> bcEp = PeriodicBC(E, PeriodicBoundary())
>
> u2_expr = Expression(("1 + x[0]*(1-x[0])","x[1]*(1-x[1])")) # curl-free
> u2 = Function(interpolate(u2_expr,S))
>
> v = TrialFunction(E)
> w = TestFunction(E)
> om2 = Function(E)
> solve(inner(v,w)*dx==-inner(as_vector([w.dx(1),-w.dx(0)]),u2)*dx,om2,bcEp)
> print om2.vector().array()
>
> As of the current dev version, the node numbering for P2 elements in UnitSquare(1,1) appears to be:
> 8 7 4
> 6 3 1
> 5 0 2
>
> The difference between the values at 1 and 6 is 0, but the difference
> between the 0-7 pair is O(10^-16). None of the corner pairs match. The
> 5-8 pair has a difference of O(10^-16), while all the other pairs have
> O(1) differences.
>
> In the UnitSquare(2,2) case, the numbering seems to be
> 24 23 19 12 11
> 22 21 18 10 08
> 20 17 16 06 03
> 15 13 07 05 01
> 14 09 04 00 02
>
> The 0-12 and 4-19 pairs weren't numerically identical, with errors at
> the 10^-16 / 10^-17 scale. Again, none of the corner pairs match, while
> the horizontal pairs 8-22, 3-20, 1-15 have differences that are
> numerically 0.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Ah, thank you. Do you have an idea when these changes will be merged into the trunk?

I previously used multiple PeriodicBC objects in the stable version, and I noticed it complained when I switched to the dev version.

Thanks Mikael Mortensen, that solved my question.

Den Nov 15, 2012 kl. 3:50 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214267
>
> Andrew McRae posted a new comment:
> Ah, thank you. Do you have an idea when these changes will be merged
> into the trunk?

Hopefully sooner than later. Meanwhile you can do like this:

from dolfin import *
set_log_level(ERROR)
mesh = UnitSquare(1, 1)

class PeriodicBoundaryX(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.)

    def map(self, x, y):
        y[0] = x[0] - 1
        y[1] = x[1]

class PeriodicBoundaryY(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.)

    def map(self, x, y):
        y[0] = x[0]
        y[1] = x[1] - 1.

left = PeriodicBoundaryX()
bottom = PeriodicBoundaryY()

E = FunctionSpace(mesh, 'CG', 2)
S = FunctionSpace(mesh, 'BDM', 1)

bcEp = [PeriodicBC(E, left), PeriodicBC(E, bottom)]

u2_expr = Expression(("1 + x[0]*(1-x[0])","x[1]*(1-x[1])")) # curl-free
u2 = Function(interpolate(u2_expr,S))

v = TrialFunction(E)
w = TestFunction(E)
om2 = Function(E)
#solve(inner(v,w)*dx==-inner(as_vector([w.dx(1),-w.dx(0)]),u2)*dx,om2)
A = assemble(inner(v,w)*dx)
b = assemble(-inner(as_vector([w.dx(1),-w.dx(0)]),u2)*dx)
[bc.apply(A, b) for bc in bcEp]
solve(A, om2.vector(), b)

print om2.vector().array()

Mikael

>
> I previously used multiple PeriodicBC objects in the stable version, and
> I noticed it complained when I switched to the dev version.

>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

This is in the stable, or in your branch? (Or both?)

This should work with the latest in trunk, not the stable one.

Mikael

Den Nov 15, 2012 kl. 5:01 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214267
>
> Andrew McRae posted a new comment:
> This is in the stable, or in your branch? (Or both?)
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Amazing, thank you!

(I didn't notice that you've split up the solve, so at first I got the same 'cannot have multiple periodic BC' error)

Hi Mikael,

Now that your work is in trunk, what is the correct way to do doubly-periodic?

Hi,

Try this for a doubly periodic domain for a rectangle mesh [-pi, pi]x[-pi, pi]

class PeriodicDomain(SubDomain):
    """Doubly periodic for rectangle [-pi, pi] x [-pi, pi]"""
    def inside(x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)
        return bool((near(x[0], -pi) or near(x[1], -pi)) and
                (not ((near(x[0], -pi) and near(x[1], pi)) or
                    (near(x[0], pi) and near(x[1], -pi)))) and on_boundary)

    def periodic_map(x, y):
        if near(x[0], pi) and near(x[1], pi):
            y[0] = x[0] - 2.0*pi
            y[1] = x[1] - 2.0*pi
        elif near(x[0], pi):
            y[0] = x[0] - 2.0*pi
            y[1] = x[1]
        elif near(x[1], pi):
            y[0] = x[0]
            y[1] = x[1] - 2.0*pi
        else:
            y[0] = 1e8
            y[1] = 1e8

Use it like in the periodic demo.

Best regards

Mikael

Den Mar 7, 2013 kl. 1:21 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214267
>
> Status: Solved => Open
>
> Andrew McRae is still having a problem:
> Hi Mikael,
>
> Now that your work is in trunk, what is the correct way to do doubly-
> periodic?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Thanks. I have it working for the CG space, but if I try to use the BDM space, I get

Traceback (most recent call last):
  File "solntest.py", line 9, in <module>
    S = FunctionSpace(mesh, 'BDM', 1, constrained_domain=pbc)
  File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 390, in __init__
    FunctionSpaceBase.__init__(self, mesh, element, constrained_domain)
  File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 89, in __init__
    dolfin_dofmap = cpp.DofMap(ufc_dofmap, mesh, constrained_domain)
  File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/cpp/fem.py", line 668, in __init__
    _fem.DofMap_swiginit(self,_fem.new_DofMap(*args))
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** https://answers.launchpad.net/dolfin
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to complete call to function build_ufc().
*** Reason: Assertion slave_master_entities->find(0) != slave_master_entities->end() failed.
*** Where: This error was encountered inside /home/atm112/fenics/src/dolfin/current/dolfin/fem/DofMapBuilder.cpp (line 426).
*** Process: 0
*** -------------------------------------------------------------------------

Is the new PeriodicBC not yet implemented for some elements?

Hi,

Could you give me a copy of the script you're running? I have no experience with BDM, but I know DG elements do not work with periodic yet. Did BDM use to work with the old PeriodicBC?

Mikael

Den Mar 7, 2013 kl. 2:20 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214267
>
> Andrew McRae posted a new comment:
> Thanks. I have it working for the CG space, but if I try to use the BDM
> space, I get
>
> Traceback (most recent call last):
> File "solntest.py", line 9, in <module>
> S = FunctionSpace(mesh, 'BDM', 1, constrained_domain=pbc)
> File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 390, in __init__
> FunctionSpaceBase.__init__(self, mesh, element, constrained_domain)
> File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 89, in __init__
> dolfin_dofmap = cpp.DofMap(ufc_dofmap, mesh, constrained_domain)
> File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/cpp/fem.py", line 668, in __init__
> _fem.DofMap_swiginit(self,_fem.new_DofMap(*args))
> RuntimeError:
>
> *** -------------------------------------------------------------------------
> *** DOLFIN encountered an error. If you are not able to resolve this issue
> *** using the information listed below, you can ask for help at
> ***
> *** https://answers.launchpad.net/dolfin
> ***
> *** Remember to include the error message listed below and, if possible,
> *** include a *minimal* running example to reproduce the error.
> ***
> *** -------------------------------------------------------------------------
> *** Error: Unable to complete call to function build_ufc().
> *** Reason: Assertion slave_master_entities->find(0) != slave_master_entities->end() failed.
> *** Where: This error was encountered inside /home/atm112/fenics/src/dolfin/current/dolfin/fem/DofMapBuilder.cpp (line 426).
> *** Process: 0
> *** -------------------------------------------------------------------------
>
> Is the new PeriodicBC not yet implemented for some elements?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Yes, BDM worked with the old PeriodicBC.

from dolfin import *
mesh = UnitSquareMesh(5, 5)
class PeriodicBoundaryX(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0)
    def map(self, x, y):
        y[0] = x[0] - 1
        y[1] = x[1]
pbc = PeriodicBoundaryX()
E = FunctionSpace(mesh, 'CG', 2, constrained_domain=pbc)
S = FunctionSpace(mesh, 'BDM', 1, constrained_domain=pbc)
#S = VectorFunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)

Ok, thanks. Looks like there's an assertion in DofMapBuilder that shouldn't be there. At least when I remove it it seems to work.

I'll make a patch as soon as I know for sure.

Mikael

Den Mar 7, 2013 kl. 3:21 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214267
>
> Status: Answered => Open
>
> Andrew McRae is still having a problem:
> Yes, BDM worked with the old PeriodicBC.
>
> from dolfin import *
> mesh = UnitSquareMesh(5, 5)
> class PeriodicBoundaryX(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and near(x[0], 0)
> def map(self, x, y):
> y[0] = x[0] - 1
> y[1] = x[1]
> pbc = PeriodicBoundaryX()
> E = FunctionSpace(mesh, 'CG', 2, constrained_domain=pbc)
> S = FunctionSpace(mesh, 'BDM', 1, constrained_domain=pbc)
> #S = VectorFunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Thanks.

By the way, is there any reason you included an 'else' clause in the code you suggested?

    def periodic_map(x, y):
        if near(x[0], pi) and near(x[1], pi):
            y[0] = x[0] - 2.0*pi
            y[1] = x[1] - 2.0*pi
        elif near(x[0], pi):
            y[0] = x[0] - 2.0*pi
            y[1] = x[1]
        elif near(x[1], pi):
            y[0] = x[0]
            y[1] = x[1] - 2.0*pi
        else:
            y[0] = 1e8
            y[1] = 1e8

On Sunday, 10 March 2013, Andrew McRae wrote:

> Question #214267 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214267
>
> Andrew McRae posted a new comment:
> By the way, is there any reason you included an 'else' clause in the
> code you suggested?
>
>
Yes, to make sure all points off boundaries are mapped off the slave
boundary. Perhaps we should try to get rid of it though? With the algorithm
used to find matching entities the map is called for entities on and off
the boundaries. So when none of the first three if statements are true y
will simply not change its value. If it happens that y was on the
slave boundary in the previous check then inside(y) will still be true. A
remedy I see now could be to set default y value in C++ before the map is
called. Then the last else should not be needed. Will look into it.

Mikael

> def periodic_map(x, y):
> if near(x[0], pi) and near(x[1], pi):
> y[0] = x[0] - 2.0*pi
> y[1] = x[1] - 2.0*pi
> elif near(x[0], pi):
> y[0] = x[0] - 2.0*pi
> y[1] = x[1]
> elif near(x[1], pi):
> y[0] = x[0]
> y[1] = x[1] - 2.0*pi
> else:
> y[0] = 1e8
> y[1] = 1e8
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

It seems to work no matter what value it's set to (including 0, 0.5, 1... which I'd otherwise expect to cause issues!). But if it's removed altogether then things break, as I just found out...