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.
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++.
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.
Thank you once again for your response. I sincerely appreciate your clever idea, which helps us a lot.
In addition to the follow-up questions mentioned in my earlier reply, I would like to ask: what potential optimizations were you referring to in your previous message?
I understand that one possible drawback of this approach is that it may lead to a larger system than expected. One idea I have considered is rewriting the system using blocks, similar to Schur complements in mixed FEM. Do you believe that further optimizations could be achieved in terms of the solver and preconditioning aspects?
Thank you in advance for any insights or suggestions you might have regarding these follow-up questions.
we typically prefer HDG instead of DG (see i-tutorials 2.7, 2.8, 3.7).
There we introduce the meanvalue on edges as independent variables, which allows for element-local assembling.
Here is an example where we build an LDG method via HDG.
A second optimization one should do is to eliminate one of the two gradient variables (already on paper).
