DG Quad elements do not converge


145
views
0
8 weeks ago by
Nate  
I'm trying to solve DG problems using the new DGQuad elements, however I can't seem to get them to converge to a known solution.

I've constructed a basic manufactured solution test for the linear Poisson equation, editing the demo code slightly.

I also assemble the cell volumes into a DG0 function since CellSize and CellVolume are currently not working with DGQuad elements.

In the following code, turn off/on the UnitSquareMesh and UnitQuadMesh to demonstrate the effect.

If anyone can spot something stupid that I've done, that would be wonderful.

from dolfin import *
import numpy as np

parameters["ghost_mode"] = "shared_facet"

run_count = 0
ele_ns = [4, 8, 16, 32, 64]
errorl2 = np.zeros(len(ele_ns))
errorh1 = np.zeros(len(ele_ns))
hsizes = np.zeros(len(ele_ns))


for ele_n in ele_ns:
  # Create mesh and define function space
  # mesh = UnitQuadMesh.create(ele_n, ele_n)
  mesh = UnitSquareMesh(ele_n, ele_n)
  V = FunctionSpace(mesh, 'DG', 1)

  # Define test and trial functions
  u = TrialFunction(V)
  v = TestFunction(V)

  # Assemble cell volumes into a DG0 function since quad elements do not support
  # CellSize() nor CellVolume()
  DG0 = FunctionSpace(mesh, "DG", 0)
  cell_volume = Function(DG0)
  v0 = TestFunction(DG0)
  cell_volume.vector()[:] = assemble(v0*dx)

  # # Define normal vector and mesh size
  n = FacetNormal(mesh)
  h = cell_volume/FacetArea(mesh)

  h_avg = (h('+') + h('-'))/2

  # Define the source term f, Dirichlet term u0 and Neumann term g
  f = Expression('2*pi*pi*sin(pi*x[0])*sin(pi*x[1])', element=V.ufl_element())
  u0 = Expression("sin(pi*x[0])*sin(pi*x[1])", element=V.ufl_element())

  # Mark facets of the mesh
  boundaries = FacetFunction('size_t', mesh, 0)
  AutoSubDomain(lambda x, on: on).mark(boundaries, 1)

  # Define outer surface measure aware of Dirichlet and Neumann boundaries
  ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

  # Define parameters
  alpha = 4.0
  gamma = 8.0

  # Define variational problem
  a = dot(grad(v), grad(u))*dx \
     - dot(avg(grad(v)), jump(u, n))*dS \
     - dot(jump(v, n), avg(grad(u)))*dS \
     + alpha/h_avg*dot(jump(v, n), jump(u, n))*dS \
     - dot(grad(v), u*n)*ds(1) \
     - dot(v*n, grad(u))*ds(1) \
     + (gamma/h)*v*u*ds(1)
  L = v*f*dx - u0*dot(grad(v), n)*ds(1) + (gamma/h)*u0*v*ds(1)

  # Compute solution
  u = Function(V)
  solve(a == L, u)

  # Assemble l2 error
  errorl2[run_count] = assemble((u0 - u)**2*dx)**0.5
  hsizes[run_count] = mesh.hmax()
  run_count += 1

# Compute convergence rates. L2 SIPG convergence rate should be ~ p+1
if dolfin.MPI.rank(mesh.mpi_comm()) == 0:
    print(np.log(errorl2[0:-1]/errorl2[1:])/np.log(hsizes[0:-1]/hsizes[1:]))​
Community: FEniCS Project
1
1
This is currently being discussed at https://bitbucket.org/fenics-project/ffc/pull-requests/91
written 5 weeks ago by Jan Blechta  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »