periodic boundary conditions

[gdmcbain]

Any ideas on periodic boundary conditions? There's a brief mention in the JOSS paper.

I'm getting increasingly interested in problems of propagation, for which these often appear in textbook examples, so I'll soon look into throwing something together.

[kinnala]

I think I used to have something earlier, let me check out the Python 2.7 version of the library.

[kinnala]

So from my archives I found something like this:

    def periodify_matrix(Z):
        Zxx=Z[Nelse].T[Nelse].T
        Zx1=Z[Nelse].T[Nleft].T
        Z1x=Zx1.T
        Zx2=Z[Nelse].T[Nright].T
        Z2x=Zx2.T
        Z11=Z[Nleft].T[Nleft].T
        Z12=Z[Nleft].T[Nright].T
        Z21=Z12.T
        Z22=Z[Nright].T[Nright].T

        Zp1=spsp.hstack((Zxx,Zx1+Zx2))
        Zp2=spsp.hstack((Z1x+Z2x,Z11+Z12+Z21+Z22))
        return spsp.vstack((Zp1,Zp2)).tocsr()

    def periodify_vector(y):
        return np.hstack((y[Nelse],y[Nleft]+y[Nright]))

    def unpack_periodic(y):
        Y=np.zeros(f.shape[0])
        Y[Nelse]=y[I1]
        Y[Nleft]=y[I2]
        Y[Nright]=y[I2]
        return Y

[kinnala]

Now while this is one way to do this, I think we could have some tools nowadays that help generalizing to arbitrary elements.

[kinnala]

Using the new Mesh class there is a fancy way of dealing with periodic meshes. I had not tried this before and I found out some rough edges, but it seems to be working. The idea is to describe the mesh not as ElementTriP1 field but using the corresponding discontinuous element:

# Should we define this in skfem.element?
class ElementTriP1DG(ElementTriP1):
    nodal_dofs = 0
    interior_dofs = 3

Now we can create a new mesh class which uses the above ElementTriP1DG as the underlying field:

from dataclasses import dataclass
from typing import Type

@dataclass
class MeshTriDG(MeshTri2):
    elem: Type[Element] = ElementTriP1DG

This allows doing magical things. Let us create the square mesh

m = MeshTri().refined()

and modify it so that the nodes on the left and on the right boundaries use the same topological indexing:

t = m.t.copy()
remap = np.arange(9, dtype=np.int64)
remap[3] = 2
remap[7] = 5
remap[1] = 0
print(t)
# [[0 1 1 2 2 3 4 6]
# [4 6 4 6 5 7 5 7]
# [5 7 6 8 6 8 6 8]]
t = remap[t]
print(t)
# [[0 0 0 2 2 2 4 6]
# [4 6 4 6 5 5 5 5]
# [5 5 6 8 6 8 6 8]]

(It took me a while to figure out what indices should be mapped but the idea is that 2 is at the left boundary and it should match 3 which is at the right boundary.)

Finally we initialize this "DG-Mesh":

M = MeshTriDG(
    m.p[:, m.t.T.flatten()],
    t,
)

And voila!

from skfem.models import laplace, unit_load

elem = ElementTriP1DG()
# rough edge: Mesh.bndelem causes Exception cause of new Element class
# therefore explicit mapping is needed
mapping = MappingIsoparametric(M, elem)
basis = CellBasis(M, ElementTriP1(), mapping)
A = laplace.assemble(basis)
f = unit_load.assemble(basis)
I = np.array([5, 6], dtype=np.int64)  # only two nodes remain in the interior
x = solve(*condense(A, f, I=I))
plot(basis, x, shading='gouraud')

[kinnala]

PS. It's mindblowing that the above worked with no modifications to the source code. Maybe we should introduce ElementTriP1DG and MeshTri1DG so that they are ready to be used?

[kinnala]

This could help:

def periodify(m, ix1, ix2, sort):
    assert len(ix1) == len(ix2)
    ix1 = ix1[sort(m.p[:, ix1])]
    ix2 = ix2[sort(m.p[:, ix2])]
    remap = np.arange(m.t.max() + 1, dtype=np.int64)
    remap[ix1] = ix2
    return remap[m.t]
    

ix1 = m.nodes_satisfying(lambda x: x[0] == 1)
ix2 = m.nodes_satisfying(lambda x: x[0] == 0)

def sort(p):
    return np.argsort(p[0])

t = periodify(m, ix1, ix2, sort)
t

[kinnala]

In PR #692 this can be done:

from skfem import *
from skfem.visuals.matplotlib import plot, draw, show()
import numpy as np

m = MeshTri().refined()
draw(m, node_numbering=True)
show()

M = MeshTri1DG.from_mesh(m)  # duplicate nodes
Mp = M.periodic(  # make periodic
    m.nodes_satisfying(lambda x: x[0] == 1),  # nodes on this set are replaced by ...
    m.nodes_satisfying(lambda x: x[0] == 0),  # ... nodes on this set
    lambda p: np.argsort(p[1]),  # couldn't yet figure if this can be automated somehow
)

from skfem.models import laplace, unit_load

basis = Basis(Mp, ElementTriP1())
A = laplace.assemble(basis)
f = unit_load.assemble(basis)
I = np.array([5, 6], dtype=np.int64)  # look up 5 and 6 from the original mesh figure
x = solve(*condense(A, f, I=I))
plot(basis, x, shading='gouraud')

[kinnala]

Above is WIP. Any feedback on, e.g., MeshTri1DG.periodic(...) is appreciated. It's a bit clumsy to use right now.

[kinnala]

Seems that {interior,boundary}_{nodes,facets} etc. do not work properly with MeshTri1DG. I need to see if anything can be done about it.

[gdmcbain]

Above is WIP. Any feedback on, e.g., MeshTri1DG.periodic(...) is appreciated.

GitHub seems to have had a glitch. I left a comment yesterday and there was a fireworks click from Bhavesh on the example, but don't see either now. Anyway my comment wasn't much more substantive than fireworks ("wow...mind blowing, ...magical, ...", &c.) but this great work is not being ignored. I'd assume I'd simply forgotten to click comment before logging off but I wouldn't think Bhavesh would have withdrawn the fireworks.

I'll take this for a spin in the morning, in particular looking for the interior/boundary nodes/facets and making substantive the speculations from my lost comment.

A minor usability idea: the sort argument might be made optional, defaulting to assuming that the user has provided the two lots of nodes in corresponding order. I think some mesh generators might provide this info; for my part though I mostly emvisage periodicity being used with meshes generated internally by scikit-fem, regular structured meshes on segment, rectangle, or cuboid.

[kinnala]

I'll take this for a spin in the morning, in particular looking for the interior/boundary nodes/facets and making substantive the speculations from my lost comment.

Don't spend too much time on it if there are any issues, it's still untested and I'm having some issues finding the correct nodes so there could be a bug lurking somewhere.

[kinnala]

Here is what I have now:

from skfem import *
from skfem.visuals.matplotlib import plot, draw
import numpy as np

m = MeshTri.init_tensor(np.linspace(0, 1, 10), np.linspace(0, 1, 10))
draw(m, node_numbering=True)

M = MeshTri1DG.from_mesh(m)
Mp = M.periodic(
    m.nodes_satisfying(lambda x: (x[0] == 1)),
    m.nodes_satisfying(lambda x: (x[0] == 0)),
)

I = m.nodes_satisfying(lambda x: (x[1] < 1 - 1e-5) * (x[0] < 1 - 1e-5) * (x[1] > 1e-5))

from skfem.models import laplace, unit_load

basis = CellBasis(Mp, ElementTriP1())
A = laplace.assemble(basis)
f = unit_load.assemble(basis)
x = solve(*condense(A, f, I=I))
plot(basis, x, shading='gouraud')

User must verify now that the parameters given to periodic are sorted properly. For this particular mesh nodes_satisfying happens to return them in correct order.

[kinnala]

Maybe there should be a new method vertices_satisfying that refers to the vertices of the underlying topology. One issue above is that we must use nodes_satisfying from the original mesh since MeshTri1DG, by design, duplicates each vertex as many times as there are elements using it. Thus, M.nodes_satisfying returns duplicates.

[kinnala]

Now basis.get_dofs() seems to work:

from skfem.models import laplace, unit_load

basis = CellBasis(Mp, ElementTriP1())
A = laplace.assemble(basis)
f = unit_load.assemble(basis)
D = basis.get_dofs()
x = solve(*condense(A, f, D=D))
plot(basis, x, shading='gouraud')

Didn't change a thing, maybe I had some typo before.

[bhaveshshrimali]

GitHub seems to have had a glitch. I left a comment yesterday and there was a fireworks click from Bhavesh on the example, but don't see either now. Anyway my comment wasn't much more substantive than fireworks ("wow...mind blowing, ...magical, ...", &c.) but this great work is not being ignored. I'd assume I'd simply forgotten to click comment before logging off but I wouldn't think Bhavesh would have withdrawn the fireworks.

I swear, I never withdrew it. GitHub indeed seems to have had a glitch. There goes the fireworks again. I can't express how excited I am at this functionality.

I'll try to give this a go sometime in the coming weeks. I am traveling Monday and busy packing today and tomorrow. And the semester starts at UIUC from Monday. But my aim is to first reproduce this example from it's corresponding implementation in FEniCS. This has been on my wishlist ever since #438 .

Many thanks @kinnala

[gdmcbain]

Oops, sorry, no, I see now that the first example had been followed by a second.

[bhaveshshrimali]

And I didn't notice either.. Anyways, both are fireworks worthy :) I'll post here when I take a crack at reproducing the above linked example.

[gdmcbain]

I found something that doesn't work: .point_source. I was thinking that the Poisson example above would better show the periodicity if the solution varied in the periodic variable and so replaced the uniform source by a point-source as in

scikit-fem/docs/examples/ex38.py

Line 46 in e4a4b2b

b = basis.point_source(source)

as in

m = MeshTri.init_tensor(np.linspace(0, 1, 10), np.linspace(0, 1, 10))

M = MeshTri1DG.from_mesh(m)
Mp = M.periodic(
    m.nodes_satisfying(lambda x: (x[0] == 1)),
    m.nodes_satisfying(lambda x: (x[0] == 0)),
)

basis = CellBasis(Mp, ElementTriP1())
f = basis.point_source(np.array([0.3, 0.2])

but the last raises

Traceback (most recent call last):
  File "/home/gdmcbain/src/scikit-fem/mwe.py", line 13, in <module>
    f = basis.point_source(np.array([0.3, 0.2]))
  File "/home/gdmcbain/src/scikit-fem/skfem/assembly/basis/cell_basis.py", line 176, in point_source
    return self.probes(x[:, None]).toarray()[0]
  File "/home/gdmcbain/src/scikit-fem/skfem/assembly/basis/cell_basis.py", line 146, in probes
    cells = self.mesh.element_finder(mapping=self.mapping)(*x)
  File "/home/gdmcbain/src/scikit-fem/skfem/mesh/mesh_tri_1.py", line 344, in finder
    X = mapping.invF(np.array([x, y])[:, None], ix)
  File "/home/gdmcbain/src/scikit-fem/skfem/mapping/mapping_isoparametric.py", line 145, in invF
    raise Exception(("Newton iteration didn't converge "
Exception: Newton iteration didn't converge up to TOL=1e-12

More directly, the same is raised with

Mp.element_finder()(*np.array([[0.3], [0.2]]))

[kinnala]

I made a change, removed periodic and now the Poisson example looks like this:

from skfem import *
from skfem.visuals.matplotlib import plot, draw
import numpy as np

m = MeshTri.init_tensor(np.linspace(0, 1, 10),
                        np.linspace(0, 1, 10))
remap = np.arange(m.nvertices, dtype=np.int64)
remap[m.nodes_satisfying(lambda x: x[0] == 1)] = m.nodes_satisfying(lambda x: x[0] == 0)
Mp = MeshTri1DG.from_mesh(m, remap[m.t])

from skfem.models import laplace, unit_load

basis = CellBasis(Mp, ElementTriP1())
A = laplace.assemble(basis)
f = unit_load.assemble(basis)
D = basis.get_dofs()
x = solve(*condense(A, f, D=D))
plot(basis, x, shading='gouraud')

[kinnala]

I found something that doesn't work: .point_source. I was thinking that the Poisson example above would better show the periodicity if the solution varied in the periodic variable and so replaced the uniform source by a point-source as in

scikit-fem/docs/examples/ex38.py

Line 46 in e4a4b2b

b = basis.point_source(source)

as in

m = MeshTri.init_tensor(np.linspace(0, 1, 10), np.linspace(0, 1, 10))

M = MeshTri1DG.from_mesh(m)
Mp = M.periodic(
    m.nodes_satisfying(lambda x: (x[0] == 1)),
    m.nodes_satisfying(lambda x: (x[0] == 0)),
)

basis = CellBasis(Mp, ElementTriP1())
f = basis.point_source(np.array([0.3, 0.2])

but the last raises

Traceback (most recent call last):
  File "/home/gdmcbain/src/scikit-fem/mwe.py", line 13, in <module>
    f = basis.point_source(np.array([0.3, 0.2]))
  File "/home/gdmcbain/src/scikit-fem/skfem/assembly/basis/cell_basis.py", line 176, in point_source
    return self.probes(x[:, None]).toarray()[0]
  File "/home/gdmcbain/src/scikit-fem/skfem/assembly/basis/cell_basis.py", line 146, in probes
    cells = self.mesh.element_finder(mapping=self.mapping)(*x)
  File "/home/gdmcbain/src/scikit-fem/skfem/mesh/mesh_tri_1.py", line 344, in finder
    X = mapping.invF(np.array([x, y])[:, None], ix)
  File "/home/gdmcbain/src/scikit-fem/skfem/mapping/mapping_isoparametric.py", line 145, in invF
    raise Exception(("Newton iteration didn't converge "
Exception: Newton iteration didn't converge up to TOL=1e-12

More directly, the same is raised with

Mp.element_finder()(*np.array([[0.3], [0.2]]))

Thanks! I think many of the methods inherited from MeshTri1 will not work. However, I needed _uniform (for plotting) and since that's quite lengthy I decided to inherit from MeshTri1 and not Mesh2D.

[kinnala]

Ok, there is an obvious "mistake" in MeshTri1.element_finder which probably also affects MeshTri2. I'll look if there is a quick fix or otherwise will disable those nonverified features.

[gdmcbain]

O. K., ta.

On a more positive note, I have replaced the point-source with a tight Gaussian and added uniform advection to the right and produced an example that does show off the periodicity: the plume from the local source is advected out the right and reappears through the left.

@BilinearForm
def advection(u, v, _):
    return v * u.grad[0]


@LinearForm
def source(v, w):
    x = w.x - 0.5
    return v * np.exp(-1e3 * (x[0]**2 + x[1]**2))


m = MeshTri.init_tensor(np.linspace(0, 1), np.linspace(0, 1))
draw(m, node_numbering=True)

M = MeshTri1DG.from_mesh(m)
Mp = M.periodic(
    m.nodes_satisfying(lambda x: (x[0] == 1)),
    m.nodes_satisfying(lambda x: (x[0] == 0)),
)

I = m.nodes_satisfying(lambda x: (x[1] < 1 - 1e-5) * (x[0] < 1 - 1e-5) * (x[1] > 1e-5))

peclet = 1e2

basis = CellBasis(Mp, ElementTriP1())
A = laplace.assemble(basis) + peclet * advection.assemble(basis)
f = source.assemble(basis)

D = basis.get_dofs()
x = solve(*condense(A, f, D=D))
ax = plot(basis, x, shading='gouraud')

[kinnala]

Cool! Today I'll work on some tests and try to have this implemented for all other mesh types as well.

[kinnala]

Took me too long to figure this out but:

1. Seems that the "eliminated" nodes should have the highest indices, probably so that there won't be gaps in the indexing. (This requires reindexing of nodes, already have a POC for that.)
2. Sometimes get_dofs is not working; seems to be unrelated to 1. Trying to figure it out.

[kinnala]

from skfem import *

m = MeshQuad1().refined(2)

m = m.with_boundaries({
    'left': lambda x: x[0] == 0,
    'right': lambda x: x[0] == 1,
})
mp = MeshQuad1DG.from_mesh(m, m.periodic_connectivity(('right', 'left')))

basis = Basis(mp, ElementQuad1())
A = laplace.assemble(basis)
f = unit_load.assemble(basis)
# D = basis.get_dofs()
# TODO get_dofs not working because f2t does not always work
I = m.nodes_satisfying(lambda x: (x[0] < 1) * (x[1] < 1) * (x[1] > 0))
x = solve(*condense(A, f, I=I))

from skfem.visuals.matplotlib import plot, draw

ax = draw(basis)
plot(basis, x, shading='gouraud', ax=ax)

[kinnala]

Still trying to figure out why get_dofs is sometimes failing.

[kinnala]

get_dofs fails because f2t is sometimes wrong. Observe:

m = MeshQuad1().refined(1).with_boundaries({
    'left': lambda x: x[0] == 0,
    'right': lambda x: x[0] == 1,
})
mp = MeshQuad1DG.from_mesh(m, m.periodic_connectivity(('right', 'left')))

print(m.f2t)
print(mp.f2t)

which returns

[[ 0  0  1  1  2  2  3  3  0  3  2  3]
 [-1 -1 -1 -1 -1 -1 -1 -1  1  0  1  2]]
[[0 1 2 2 0 2 3]
 [1 0 3 3 1 1 2]]

So mp has no boundary facets? Obviously wrong. However, if you refine once more it suddenly has boundary facets.

Have not yet found out what causes this. For some meshes it works, for some meshes it doesn't.

[kinnala]

Now things like this should work:

from skfem import *
import numpy as np
from skfem.models import laplace, unit_load, mass
from skfem.visuals.matplotlib import draw, plot

m = MeshTri().refined(2)
mp = MeshTri1DG.periodic(
    m,
    m.nodes_satisfying(lambda x: (x[0] == 1)),
    m.nodes_satisfying(lambda x: (x[0] == 0)),
)

e = ElementTriP1()
basis = Basis(mp, e)
A = laplace.assemble(basis)
f = unit_load.assemble(basis)
x = solve(*condense(A, f, D=basis.get_dofs().flatten()))

ax = draw(basis)
plot(basis, x, shading='gouraud', ax=ax)

This constructor, Mesh.periodic, does not allow the creation of more generic periodic meshes (e.g. in two directions). However, that should still work if the mesh is created manually. I'll try that out shortly.

[kinnala]

Note: higher order elements work:

from skfem import *
import numpy as np
from skfem.models import laplace, unit_load, mass
from skfem.visuals.matplotlib import draw, plot

m = MeshQuad().refined(2)
mp = MeshQuad1DG.periodic(
    m,
    m.nodes_satisfying(lambda x: (x[0] == 1)),
    m.nodes_satisfying(lambda x: (x[0] == 0)),
)

e = ElementQuadP(6)
basis = Basis(mp, e)
A = laplace.assemble(basis)
f = unit_load.assemble(basis)
x = solve(*condense(A, f, D=basis.get_dofs().flatten()))

plot(basis, x, shading='gouraud', nrefs=3)

[kinnala]

Setting this up was a bit tedious but it seems to work. Note that I removed the sorting inside periodic so the user must verify that the arguments are sorted:

from skfem import *
import numpy as np
from skfem.models import laplace, unit_load, mass
from skfem.visuals.matplotlib import draw, plot

m = MeshQuad().refined(5)

ix1 = m.nodes_satisfying(lambda x: (x[0] == 1) * (x[1] < 1) * (x[1] > 0))
ix2 = m.nodes_satisfying(lambda x: (x[0] == 0) * (x[1] < 1) * (x[1] > 0))

ix1 = np.concatenate((
    ix1,
    m.nodes_satisfying(lambda x: (x[1] == 0) * (x[0] < 1) * (x[0] > 0)),
))
ix2 = np.concatenate((
    ix2,
    m.nodes_satisfying(lambda x: (x[1] == 1) * (x[0] < 1) * (x[0] > 0)),
))

ix1 = np.concatenate((
    ix1,
    m.nodes_satisfying(lambda x: (x[0] == 1) * (x[1] == 1)),
))

ix2 = np.concatenate((
    ix2,
    m.nodes_satisfying(lambda x: (x[0] == 0) * (x[1] == 0)),
))

ix1 = np.concatenate((
    ix1,
    m.nodes_satisfying(lambda x: (x[0] == 1) * (x[1] == 0)),
))

ix2 = np.concatenate((
    ix2,
    m.nodes_satisfying(lambda x: (x[0] == 0) * (x[1] == 1)),
))
    

mp = MeshQuad1DG.periodic(
    m,
    ix1,
    ix2,
)

e = ElementQuad1()
basis = Basis(mp, e)
A = laplace.assemble(basis)
M = mass.assemble(basis)
    
@LinearForm
def load(v, w):
    return np.sin(2. * np.pi * (w.x[0] - 0.2)) * np.sin(2. * np.pi * (w.x[1] - 0.2))

f = load.assemble(basis)
x = solve(*condense(A, f, D=basis.get_dofs().flatten()))

plot(basis, x, shading='gouraud', nrefs=3)

[kinnala]

I think we might need reverse mapping to the original mesh to do some saving etc.

[gdmcbain]

On saving, I'm not sure how important it is, but thought I'd mention it while I think of it and as there has been continued reference to the Gmsh file formats in #680 and #696:

Both Gmsh file formats, viz. MSH ‘current’ 4.1 and ‘legacy’ 2.2, contain $Periodic blocks (though not the same). It looks like meshio.gmsh does _read_periodic from both formats, but I haven't tried it. And there is _write_periodic too.

I imagine a lot of uses of periodicity (all mine that I currently envisage, and I imagine all those related to homogenization as in #438) will only involve line segment, rectangle, or cuboid domains and so not require external meshing?

[kinnala]

I suppose adding some NotImplementedError's is required.

[kinnala]

Added MeshDG.load/save with explicit NotImplementedError. Also MeshDG.element_finder raises NotImplementedError.