DG in a periodic domain

Hi,

Den Nov 29, 2012 kl. 12:11 PM skrev Stephan Kramer:

> New question #215542 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/215542
>
> Hi,
>
> I would like to solve a DG problem in a periodic domain, i.e. have interior facet integrals dS where the facets on the periodic boundary are simply treated as internal, and stuff like jump() and avg() would work as well. I understand that this is not currently supported in Dolfin where periodic boundary conditions are treated as strong dirichlet bcs. I saw the branch lp:~mikael-mortensen/dolfin/periodic_to_dofmap/ and was hoping that would implement it, as it seems to work at the mesh level - so it could in principle fix the topology to be periodic so that all derived function spaces are automatically periodic.

You are correct and I made a note of this earlier when I suggested adding the periodicity directly to the mesh.

>
> So my question really is, is this indeed how the new periodic boundary conditions on that branch work, and would periodic DG interior facet integrals indeed be supported? I tried out the branch and it didn't work as expected. See code below, using the new add_periodic_direction() method, for a simple scalar advection equation, on a domain periodic in the x-direction, advecting a scalar field that only varies in the y-direction, and thus should remain unchanged (but it doesn't). Please let me know if I'm simply being impatient, or whether this is not something that's planned to be supported for the moment.

I think it should be relatively easy to implement, but this needs to be added on the mesh connectivity/topology level and currently I have focused only on the dofmaps. The facets on the periodic boundaries need to have two adjacent cells, like any other interior facet. If that is introduced I think the DG-methods will work out of the box as well.

Does anyone know how this could be entered into the mesh topology? There is already a periodic facet-to-facet map computed by the mesh class in lp:~mikael-mortensen/dolfin/periodic_to_dofmap/.

Mikael

>
> Best wishes
> Stephan Kramer
>
> from dolfin import *
> mesh = UnitSquareMesh(3,3)
> dt=0.01
>
> # Sub domain for Periodic boundary condition
> class PeriodicBoundary(SubDomain):
>
> # Left boundary is "target domain" G
> def inside(self, x, on_boundary):
> return abs(x[0]) < DOLFIN_EPS and on_boundary
>
> # Map right boundary (H) to left boundary (G)
> def map(self, x, y):
> y[0] = x[0] - 1.0
> y[1] = x[1]
>
> pbc=PeriodicBoundary()
> mesh.add_periodic_direction(pbc)
> V = FunctionSpace(mesh, "DG", 1)
> v = TestFunction(V)
> phi_new = TrialFunction(V)
> u = Expression(("1.0", "0.0")) # constant velocity
> ic = Expression("1.0") # initial solution, only varies in y
> phi = interpolate(ic, V)
>
> n = FacetNormal(mesh)
> un = (dot(u, n) + abs(dot(u, n)))/2.0
> F = v*(phi_new-phi)/dt*dx + \
> - dot(grad(v), u)*phi*dx \
> + dot(jump(v), un('+')*phi('+') - un('-')*phi('-') )*dS \
> + dot(v, un*phi)*ds
>
> for i in range(20):
> solve(lhs(F)==rhs(F), phi)
>
> plot(phi, interactive=True)
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.