Discontinuous galerkin with parameter varying across subdomains

mathematics
1

I am solving the problem

-\kappa\nabla^2u=f in \Omega

with Dirichlet boundary conditions on left and right boundaries and Neumann boundary conditions on top, bottom and middle boundaries. The domain \Omega is made of \Omega_{\mathrm{left}} and \Omega_{\mathrm{right}} and \kappa varies with subdomain such that \kappa_{\mathrm{left}}=0.1 and \kappa_{\mathrm{right}}=1.

The goal is to enforce a discontinuity in u in the middle boundary while also allowing for variation of \kappa in each subdomain. I’m using Discontinuous Galerkin with interior penalty numerical fluxes (modified).

We have, at the middle boundary
-\kappa \nabla u = D0 [[u]]
where [[u]] is jump(u, n).

I am obtaining the right jump in the middle boundary, but I am not obtaining variation of the slope in each subdomain. Can someone help debug my formulation below:

D0 = 4000
U0 = ufl.as_vector((0.66, 0.66, 0.66))

α = 100
γ = 1e5

F = inner(Κ * grad(u), grad(v)) * dx

# Add DG/IP terms
F += - inner(avg(grad(v)), jump(u, n)) * dS
F += - inner(jump(v, n), avg(grad(u))) * dS
F +=  γ / h_avg * inner(jump(v, n), jump(u, n)) * dS

# Nitsche boundary terms
F += - Κ * v * inner(grad(u), n) * ds(left)
F += - Κ * v * inner(grad(u), n) * ds(right)
F += - g * v * ds(insulated)
F += - Κ * (u - u_left) * inner(grad(v), n) * ds(left)
F += - Κ * (u - u_right) * inner(grad(v), n) * ds(right) 
F += -h / α * u * inner(grad(v), n) * ds(insulated)
F += α / h * (u - u_left) * v * ds(left) 
F += α / h * (u - u_right) * v * ds(right)

# Middle boundary
F += + inner(avg(grad(v)), 1 / D0 * (Κ * grad(u))('-') - U0) * dS(middle)
F += + inner(jump(v, n), avg(Κ * grad(u))) * dS(middle)
F += - γ / h_avg * inner(jump(v, n), 1 / D0 * (Κ * grad(u))('+') - U0) * dS(middle)

# Source term
F += -f * v * dx

The rest of the MWE code is below, and includes mesh generation steps.

import os

import dolfinx
import gmsh
import matplotlib.pyplot as plt
import meshio
import numpy as np
import ufl
import warnings

from dolfinx import cpp, default_scalar_type, fem, io, mesh, nls, plot
from dolfinx.fem import petsc
from dolfinx.io import VTXWriter
from dolfinx.nls import petsc as petsc_nls
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction,
                 dot, div, dx, ds, dS, grad, inner, grad, avg, jump)

import commons, geometry, utils

warnings.simplefilter('ignore')


def create_mesh(mesh, cell_type, prune_z=False):
    """
    """
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points,
                           cells={cell_type: cells},
                           cell_data={"name_to_read": [cell_data]}
                           )
    if prune_z:
        out_mesh.prune_z_0()

    return out_mesh

if __name__ == '__main__':
    # facet labels
    left = 0
    right = 2
    middle = 3
    insulated = 4

    # subdomain labels
    left_domain = 0
    right_domain = 1

    encoding = io.XDMFFile.Encoding.HDF5
    LX = 50e-6
    LY = 250e-6
    points = [
        (0, 0, 0),
        (0.5 * LX, 0, 0),
        (LX, 0, 0),
        (LX, LY, 0),
        (0.5 * LX, LY, 0),
        (0, LY, 0),
    ]
    gpoints = []
    lines = []

    gmsh.initialize()
    gmsh.model.add('full-cell')
    # gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.1 * micron)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.5e-6)
    for idx, p in enumerate(points):
        gpoints.append(
            gmsh.model.occ.addPoint(*p)
        )
    gmsh.model.occ.synchronize()
    gmsh.model.occ.synchronize()
    for idx in range(0, len(points)-1):
        lines.append(
            gmsh.model.occ.addLine(gpoints[idx], gpoints[idx+1])
        )
    lines.append(
        gmsh.model.occ.addLine(gpoints[-1], gpoints[0])
    )
    lines.append(
        gmsh.model.occ.addLine(gpoints[1], gpoints[4])
    )

    gmsh.model.occ.synchronize()
    ltag = gmsh.model.addPhysicalGroup(1, [lines[-2]], left)
    gmsh.model.setPhysicalName(1, ltag, "left")
    rtag = gmsh.model.addPhysicalGroup(1, [lines[2]], right)
    gmsh.model.setPhysicalName(1, rtag, "right")
    evptag = gmsh.model.addPhysicalGroup(1, [lines[idx] for idx in [-1]], middle)
    gmsh.model.setPhysicalName(1, evptag, "middle")
    insultag = gmsh.model.addPhysicalGroup(1, [lines[idx] for idx in [0, 1, 3, 4]], insulated)
    gmsh.model.setPhysicalName(1, insultag, "insulated")
    gmsh.model.occ.synchronize()
    left_loop = gmsh.model.occ.addCurveLoop([lines[idx] for idx in [0, 6, 4, 5]])
    right_loop = gmsh.model.occ.addCurveLoop([lines[idx] for idx in [1, 2, 3, 6]])
    gmsh.model.occ.synchronize()
    left_phase = gmsh.model.occ.addPlaneSurface([left_loop])
    right_phase = gmsh.model.occ.addPlaneSurface([right_loop])
    gmsh.model.occ.synchronize()
    gmsh.model.addPhysicalGroup(2, [left_phase], left_domain)
    gmsh.model.addPhysicalGroup(2, [right_phase], right_domain)
    gmsh.model.occ.synchronize()

    gmsh.model.mesh.generate(2)
    gmsh.write("mesh.msh")
    gmsh.finalize()

    mesh_2d = meshio.read("mesh.msh")
    tria_mesh = create_mesh(mesh_2d, "triangle")
    meshio.write("triangles.xdmf", tria_mesh)
    line_mesh = create_mesh(mesh_2d, "line")
    meshio.write("lines.xdmf", line_mesh)

    # simulation
    comm = MPI.COMM_WORLD
    with io.XDMFFile(comm, "triangles.xdmf", "r") as infile3:
        domain = infile3.read_mesh(cpp.mesh.GhostMode.none, 'Grid')
        ct = infile3.read_meshtags(domain, name="Grid")
    tdim = domain.topology.dim
    fdim = tdim - 1
    domain.topology.create_connectivity(tdim, fdim)

    ft_imap = domain.topology.index_map(fdim)
    num_facets = ft_imap.size_local + ft_imap.num_ghosts
    indices = np.arange(0, num_facets)
    values = np.zeros(indices.shape, dtype=np.intc)  # all facets are tagged with zero

    with io.XDMFFile(comm, "lines.xdmf", "r") as infile2:
        ft = infile2.read_meshtags(domain, name="Grid")

    values[ft.indices] = ft.values
    meshtags = mesh.meshtags(domain, fdim, indices, values)
    domaintags = mesh.meshtags(domain, domain.topology.dim, ct.indices, ct.values)

    dx = ufl.Measure("dx", domain=domain, subdomain_data=domaintags)
    ds = ufl.Measure("ds", domain=domain, subdomain_data=meshtags)
    dS = ufl.Measure("dS", domain=domain, subdomain_data=meshtags)
    V = fem.FunctionSpace(domain, ("DG", 1))
    u = fem.Function(V)
    u.name = "potential"
    v = ufl.TestFunction(V)
    n = ufl.FacetNormal(domain)
    x = ufl.SpatialCoordinate(domain)

    h = ufl.CellDiameter(domain)
    h_avg = avg(h)

    f = fem.Constant(domain, PETSc.ScalarType(0))
    g = fem.Constant(domain, PETSc.ScalarType(0))
    voltage = 1000e-3
    u_left = fem.Function(V)
    with u_left.vector.localForm() as u0_loc:
        u0_loc.set(0)
    u_right = fem.Function(V)
    with u_right.vector.localForm() as u1_loc:
        u1_loc.set(voltage)

    # domain-varying parameter
    Κ = fem.Function(V)
    cells_left = ct.find(left_domain)
    cells_right = ct.find(right_domain)
    Κ.x.array[cells_left] = np.full_like(cells_left, 0.1, dtype=default_scalar_type)
    Κ.x.array[cells_right] = np.full_like(cells_right, 1, dtype=default_scalar_type)

    # constant for kinetics with -K * inner(grad(u), n) = D0 * (jump(u, n) - U0)
    D0 = 4000
    U0 = ufl.as_vector((0.66, 0.66, 0.66))

    α = 100
    γ = 1e5

    F = inner(Κ * grad(u), grad(v)) * dx

    # Add DG/IP terms
    F += - inner(avg(grad(v)), jump(u, n)) * dS
    F += - inner(jump(v, n), avg(grad(u))) * dS
    F +=  γ / h_avg * inner(jump(v, n), jump(u, n)) * dS

    # Nitsche boundary terms
    F += - Κ * v * inner(grad(u), n) * ds(left)
    F += - Κ * v * inner(grad(u), n) * ds(right)
    F += - g * v * ds(insulated)
    F += - Κ * (u - u_left) * inner(grad(v), n) * ds(left)
    F += - Κ * (u - u_right) * inner(grad(v), n) * ds(right) 
    F += -h / α * u * inner(grad(v), n) * ds(insulated)
    F += α / h * (u - u_left) * v * ds(left) 
    F += α / h * (u - u_right) * v * ds(right)

    # Middle boundary
    F += + inner(avg(grad(v)), 1 / D0 * (Κ * grad(u))('-') - U0) * dS(middle)
    F += + inner(jump(v, n), avg(Κ * grad(u))) * dS(middle)
    F += - γ / h_avg * inner(jump(v, n), 1 / D0 * (Κ * grad(u))('+') - U0) * dS(middle)

    # Source term
    F += -f * v * dx

    # problem and solution
    problem = petsc.NonlinearProblem(F, u)
    solver = petsc_nls.NewtonSolver(comm, problem)
    solver.convergence_criterion = "residual"
    solver.maximum_iterations = 5

    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "preonly"
    opts[f"{option_prefix}pc_type"] = "lu"
    ksp.setFromOptions()
    n_iters, converged = solver.solve(u)

    with VTXWriter(comm, "potential.bp", [u], engine="BP4") as vtx:
        vtx.write(0.0)
2

Sadly I can’t run your code due to version clashes, and unfortunately there’s a lot to unpack here in your question. DG methods are eminent in their facility to weakly enforce fluxes between facets, however, the choice of these fluxes profoundly affects the numerical properties of the discretisation. I’d suggest discussing the topic with your supervisor or collaborators, or conducting a literature review following from the famous paper: https://www-users.cse.umn.edu/~arnold/papers/dgerr.pdf.

I can make some comments:

  • From what you’ve written it’s not obvious to me what problem you actually want to solve. The constraint you’ve stated looks like you’re essentially demanding continuity of the flux across the interior interface. But I’m not sure what the context is and therefore it’s difficult to give advice.
  • Glancing at your code I notice you have two interior penalty parameters, alpha = 100 and gamma = 1e5. Your choice of gamma being so large makes me concerned that you’re effectively using a penalty method, not a DG method.
  • Perhaps if you refer to your problem as written in the literature I can take a look if I find time.
3

Thank you for your response.

The context of the problem is charge transfer across an interface where you have flux continuity (charge/material conservation) and a potential jump across the interface. I can solve the problem when the conductivity of the two subdomains are the same, which makes me think that perhaps I have \kappa where it should not be in my bilinear form.

I have combed through the famous Arnold paper on unified DG, and have a good idea on what properties the numerical fluxes ought to have. The problem posed in the Arnold paper is -\nabla^2u=f in \Omega subject to the boundary conditions. I realize that my problem has two concerns (1) allowing variation of \kappa in subdomains and (2) enforcing discontinuity in the field u while guaranteeing continuity of \nabla u. I can solve for (1) using Continuous Galerkin, but the final problem—of which I extracted this subproblem—is to use Discontinuous Galerkin. Separately, I can solve (2), therefore (1) as an isolated problem is what’s important. I will edit the problem title and description to capture this. In summary, what I need is to figure out how to reimplement Defining subdomains for different materials — FEniCSx tutorial in Discontinuous Galerkin.

I am adding details of my formulation below in case someone notices what I miss.

That said, I was under the impression that I could use interior penalty for problem (1) as described with the parameter \kappa carrying over into IP numerical fluxes such that:

\hat{u}_K = \langle u_h \rangle and \hat{\sigma}_K = \langle \kappa \nabla_h u_h \rangle - \alpha \langle\kappa\rangle[[u_h]]

where \langle . \rangle denotes average and [[.]] denotes jump.

I went over the motions of deriving primal formulation beginning with a restatement of my problem into a first order system.

Before plugging in my choice of numerical fluxes, I have the bilinear form:

B_h(u_h, v) = \int_{\Omega}\kappa \nabla_u u_h \cdot \nabla_h v \mathrm{dx} + \int_{\Gamma}\left( [[\kappa \hat{u} - \kappa u_h]] \cdot \langle \nabla_h v\rangle - \hat{\sigma}\cdot [[v]]\right) \mathrm{ds} + \int_{\Gamma^{\circ}}\left( \langle \kappa \hat{u} -\kappa u_h\rangle [[\nabla_h v]] - [[\hat{\sigma}]]\langle v \rangle \right) \mathrm{ds}

where \Gamma indicates all facets and \Gamma^{\circ} indicates interior facets.

The term integrated over \Gamma^{\circ} is zero due to consistency and conservativity of numerical fluxes.

Before plugging in my numerical fluxes, I have:

B_h(u_h, v) = \int_{\Omega}\kappa \nabla_u u_h \cdot \nabla_h v \mathrm{dx} + \int_{\Gamma}\left( [[\kappa \hat{u} - \kappa u_h]] \cdot \langle \nabla_h v\rangle - \hat{\sigma}\cdot [[v]]\right) \mathrm{ds}

I used the above to generate my dolfinx formulation.