Implementing local discontinuous Galerkin method and assembling stiffness matrix in python interface

NGSolve
1

Hi,

This post is a follow-up to my previous post found here: Implementing Lifting Operator in DG.

I am currently working on implementing a local discontinuous Galerkin (LDG) method for the Laplace equation with homogeneous Dirichlet boundary conditions. The method and formulation are described in the following image.

LDG-4
LDG-41668×2157 534 KB

A key point of this method involves a direct use of the discrete gradient operator in the discrete formulation, which requires solving local problems (5) to compute the lifting operators.

I have written two codes using the NGSolve Python interface, which I’ve attached for reference.
LDG_new.py (12.4 KB)
LDG_poisson.py (4.7 KB)
In LDG_poisson.py, I solve the lifting operator globally, mainly using the symbolic integrator in NGSolve. Moreover, LDG_new.py solves the lifting operator more locally, utilizing various NumPy functions. While both codes are working, they are quite inefficient.

I would greatly appreciate any suggestions for improving the efficiency of the code. I am curious whether it is feasible to achieve reasonable performance using the Python interface, which is quite convenient, or if it would be more effective to transition to the base code in C++.

Thank you in advance for your suggestions!

2

Hi Syang,

you can specify the whole method including the discrete gradient via a variational formulation:

LDG.ipynb (4.1 KB)

it’s a starting point, further optimizations are possible
Joachim

3

Hi Joachim,

Many thanks. This is helpful.

I noticed that you used dS = dx(element_boundary=True). May I know what is the difference between dx(skeleton=True) and dx(element_boundary=True) in ngsolve?

Additionally, do you believe that computing the lifting operator and using it to directly assemble the stiffness matrix efficiently could be difficult within the Python interface of NGSolve? While your code is significantly faster than my previous implementation, the resulting stiffness matrix seems to be quite large.

Best