Refinement issue when use refine() in a DG formulation

mesh
1

Hi everyone,
I am trying to get a local solution from a DG implementation. To do that, I require only the inner information (’+’) in the interior facet integration (i.e. dS) instead of using the jumps and averages.
The issue is, apparently, the facet normal directions (’+’) change when I use the function refine(), but not when I do a manual uniform refinement (creating a new refined mesh on each level).

To illustrate, I wrote this short code, where I calculate the L2 norm using 1) a manual uniform refinement and 2) using the refine() function.

For the first case I have a convergent behaviour in the L2 norm
l2norm: [19.0606515, 1.6594006, 0.5111520, 0.37043388, 0.331217411, 0.316423740].

But in the second case (same code just commenting the line 37), I have the results I want only for the first two solutions (first refinement), because then the global matrix is not the same anymore.
l2norm: [19.06065, 1.6594006, inf, nan, inf, inf].

If you have any suggestions on how to fix the problem, I would appreciate it!

Thank you,
Juan

CODE:

1. import fenics
2. from dolfin import *
3. import numpy as np

4. def bh1(u,v):
5.     val = inner(kappa*grad(u), grad(v))*dx     
6.     val += -inner(u, dot(kappa*grad(v), n))*ds
7.     val += -inner(v, dot(kappa*grad(u), n))*ds
8.     val += inner((eta/h*kappa)*v, u)*ds
9.     val += dot(beta_vec, grad(u))*v*dx
10.     val += 0.5*(abs(dot(beta_vec,n))-dot(beta_vec,n))*u*v*ds
11.     return val

12. def bhfacet(u,v):   
13.     val = -inner(dot(avg(kappa*grad(u)), n('+')), jump(v))*dS 
14.     val += -inner(jump(u), dot(avg(kappa*grad(v)), n('+')))*dS
15.     val += inner((avg(eta/h)*avg(kappa))*jump(u), jump(v))*dS  
16.     val -= dot(beta_vec('+'),n('+'))*jump(u)*avg(v)*dS
17.     val += 0.5*abs(dot(beta_vec('+'),n('+')))*jump(u)*jump(v)*dS
18.     return val

19. def bhfacetmod(u,v):
20.     val = -inner(dot(kappa('+')*grad(u)('+'), n('+')), v('+'))*dS 
21.     val += -inner(u('+'), dot(kappa('+')*grad(v)('+'), n('+')))*dS
22.     val += inner((eta/h('+')*kappa('+'))*u('+'), v('+'))*dS   
23.     val -= dot(beta_vec('+'),n('+'))*u('+')*v('+')*dS
24.     val += 0.5*abs(dot(beta_vec('+'),n('+')))*u('+')*v('+')*dS
25.     return val

26. def bhdg(u,v):
27.    return bh1(u,v) + bhfacet(u,v)
28. def bhdglocal(u,v):
29.    return bh1(u,v) + bhfacetmod(u,v)

30. kappa  = Constant(1e-2) ; beta_vec= as_vector([1,]); eta = 18
31. nlist    = [40,80,160,320,640,1280] 
32. l2norm   = []
33. mesh = IntervalMesh(nlist[0], 0, 1)
34. u_exa = Expression('(exp(1/1e-2*x[0])-1)/(exp(1/1e-2)-1)-x[0]',degree=5)

35. MAX_ITER = np.size(nlist)
36. for level in range(MAX_ITER):       
37.     mesh = IntervalMesh(nlist[level], 0, 1)
38.     V = FunctionSpace(mesh, 'DG', 1)
39.     
40.     u = TestFunction(V)
41.     v = TrialFunction(V) 
42.     
43.     boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
44.     boundary_markers.set_all(0)    
45.     ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
46.     
47.     n = FacetNormal(mesh); h = CellDiameter(mesh)
48.     L = inner(v, -Constant(1))*dx
49.     
50.     solDG = Function(V)
51.     solve(bhdg(u,v) == L ,solDG )
52.     
53.     bdrhs= bhdglocal(solDG,v)
54.     
55.     uDGmod = Function(V)
56.     solve(bhdglocal(u,v) == L + bdrhs ,uDGmod)

57.     l2norm.append(assemble(inner(u_exa-uDGmod,u_exa-uDGmod)*dx))
58.     mesh = refine(mesh)
59.     
60.     
61. print('l2norm:',l2norm)
2

You cant expect the “+” restriction to be consistent, see: Integrating over an interior surface - #4 by MiroK

I guess you want to integrate over every cell, and then for each facet of the cell compute the dS integrals. This would be that you would integrate over each facet twice (one from each side). Is that what you want?

3

Thanks, Dokken, for your response.
I want to integrate over every cell and for each facet of that cell, compute only the internal “+” component. The external component (’-’) is compensated as a forcing in the RHS… So this is like I would solve the local systems.
The issue is that It works very well when I do a uniform refinement manually (creating a new mesh every single loop) because, I assume, it knows naturally which is the “+” and “-” restriction in a new mesh. But when I refine (using refine()), it is not consistent anymore after the third level of refinement.
Following the link you sent me, I assume I must indicate a restriction over all the interior surfaces every time I refine. Still, it is unclear how to indicate it because the restriction is not for a specific node but over the whole domain.
Thank you so much.

4

To integrate over every internal facet of every cell, you would have to compute both a ("+") restriction and ("-") restriction, i.e. say you have an integrand

I = inner(f, u)

you should have the form:

L = (I("+") + I("-"))*dS
5

That is clear to me. But in my case, I am not trying to solve a global DG formulation but a local one. For instance, if I have a mesh of 10 elements, it would be equivalent to solving 10 disconnected elements, and each element contains different Dirichlet boundary conditions (which I impose in the RHS). That is why I want to get only the internal component (’+’) (now ‘dS’ would represent the ‘ds’ in a global formulation).
When I create a new mesh, it is clear that fenics knows that (’+’) is the internal component (that’s why I get convergence imposing a manual refinement). But when I refine it (using refine()), it is not the internal component anymore (I just saw it in this discussion: Bitbucket). But If I can redefine what (’+’/’-’) is in a refined mesh, it would be perfect for me. If you think it is possible, I would appreciate it; thanks!

6

This is a false assumption. The plus and minus restriction is arbitrary, as discussed in the post by mirok above.
For some meshes and for a certain number of processes, you might get lucky, and you get that this restriction is what you want. As long as you do not connect the (“+”) and (“-”) restriction in the same integral (as I did above), the problem is fully disconnected (no coupling between the cells).

7

That is exactly what I want, a fully disconnected system where dS(’+’) becomes ds for each cell (and assuming, for instance, zero Dirichlet BC on each element). But it is Ok, if (’+’) and (’-’) are entirely arbitrary, it is clear that I cannot use them for my porpoise :). However, if you may suggest to me another idea to solve this disconnected system, (without having to create multiple subdomains), I would appreciate it! Thank you so much for your help

8

The point is that you need both ("+") and ("-") in your variational form to loop over each cell (and then each facet).

Say you have a 1x1 unit square mesh, then

from dolfin import *

mesh = UnitSquareMesh(1, 1)
c = Constant(1)
dS = Measure("dS", domain=mesh)
print(assemble(c("+")*dS))
print(assemble(c("-")*dS))

gives you

1.4142135623730951
1.4142135623730951

which means that this facet is only integrated once for ("+") and once for ("-")

As you want to integrate this facet twice, once for each cell connected to the facet, you need to integrate over both. This can be further illustrated by creating a DG 0 function, that is 3 in one cell, 7 in the other:

V = FunctionSpace(mesh, "DG", 0)
u = Function(V)
u.vector()[0] = 3
u.vector()[1] = 7
print(assemble(u("+")*dS))
print(assemble(u("-")*dS))
print(3*sqrt(2), 7*sqrt(2))

yielding

4.242640687119286
9.899494936611665
4.242640687119286 9.899494936611665
9

Ok thanks, Dokken, I got your point. Evidently, this is not the best way to solve this problem. I’ll assemble the matrix element by element then! thanks