MsFEM.jl

This repository contains code implementing the high-order multiscale finite element method in one and two space dimensions. The goal is to provide a more maintainable, easy-to-run codebase for Cartesian meshes. This repository will eventually be merged into the main repository (MultiscaleFEM.jl). The code uses Gridap.jl to define the Cartesian model, FESpaces and the bilinear/linear forms. Below are the instructions for running the code. The same steps can be applied to both 1D and 2D codes. First cd into the project directory (pLOD1d or pLOD2d), open the Julian prompt and instantiate the environments:

using Pkg
Pkg.instantiate()

A helper function parse_command_line() written using ArgParse.jl is provided to accept the discretization parameters from the command line. The following Julia snippet is used to extract the command line arguments:

parsed_args = parse_command_line()

n = parsed_args["fine_scale"]
N = parsed_args["coarse_scale"]
p = parsed_args["order"]
l = parsed_args["patch_radius"]
j = parsed_args["correction_level"]
ref_sol = parsed_args["reference_sol"]
out_file = parsed_args["output_file"]

Do

julia --project=. [SCRIPT] --help

to get the help environment. Once done setting up, there are two main scripts in the folders:

• poisson.jl: Solution to the Poisson problem using the p-LOD and sp-LOD methods.
• wave_equation.jl: Solution to the wave equation using p-LOD, sp-LOD and eho-LOD methods.

Explanation

pLOD vs spLOD

Let us consider the example of solving the following 1D BVP

$$ \begin{align*} &-\frac{d}{dx}\left( A(x) \frac{du}{dx} \right) = f(x), \quad x \in (0,1),\\ &u = 0, \quad x \in \{0, 1\}, \end{align*} $$

and assume that $A(x)$ oscillates at the scale $\varepsilon \ll 1$. We are interested to solve the multiscale problem using the pLOD [1] and spLOD methods [2],[3]. Run

cd pLOD1d/
julia --project=. poisson.jl -n 2048 -N 128 -p 1 -l 2
julia --project=. poisson.jl -n 2048 -N 128 -p 1 -l 7
julia --project=. poisson.jl -n 2048 -N 128 -p 1 -l 128

Along with some progress bars, we obtain the following errors:

───────────────────────────────────────────────────────────────────────────────────
  1/h  1/H  p  l                 pLOD          spLOD, HLM26           spLOD, DHM23 
───────────────────────────────────────────────────────────────────────────────────
 2048  128  1  2  0.14161089692859816  0.005022821275581407  0.0021884279093926707
 2048  128  1  7  3.916610047129179e-5  3.808523833296026e-7  2.403394931851059e-7
 2048  128  1  128  2.6309649981017746e-8  2.6309651750208505e-8  2.63096528929799e-8
──────────────────────────────────────────────────────────────────────────────────────

As proposed in [2],[3], we see that the solution obtained using the sp-LOD method is orders of magnitude better than the p-LOD method when patch radius $l=2, 7$, and they are equal (close to machine precision) when $l=\infty$. The two strategies produce slightly different solutions, showing that they generate slightly different spaces, but the errors are of similar size.

sp-LOD vs eho-LOD

Consider the 1D initial boundary value problem

$$ \begin{align*} &\frac{\partial^2 u}{\partial t^2} - \frac{\partial}{\partial x}\left( A(x) \frac{\partial u}{\partial x} \right) = f(x, t), \quad x \in (0,1), \; t > 0,\\ &u = 0, \quad x \in \{0, 1\}, \; t > 0,\\ &u(x,0) = 0, \quad x \in (0,1),\\ &\frac{\partial u}{\partial t}(x,0) = 0, \quad x \in (0,1), \end{align*} $$

and assume that $A(x)$ oscillates at the scale $\varepsilon \ll 1$. We are interested in demonstrating the effects of the additional correction bases of the eho-LOD method proposed in [4],[5]. Run the following commands:

julia --project=. wave_equation.jl -n 2048 -N 128 -p 1 -l 5 -j 0
julia --project=. wave_equation.jl -n 2048 -N 128 -p 1 -l 5 -j 1
julia --project=. wave_equation.jl -n 2048 -N 128 -p 1 -l 128 -j 0
julia --project=. wave_equation.jl -n 2048 -N 128 -p 1 -l 128 -j 1

Again, we obtain the following errors:

──────────────────────────────────────────────────────────────────────────────────────────
  1/h  1/H  p  l  j                    pLOD           spLOD, DHM23           spLOD, HLM26 
──────────────────────────────────────────────────────────────────────────────────────────
 2048  128  1  5  0  0.00028208625863618264  1.3205057333849196e-6  1.7837515480308695e-6
 2048  128  1  5  1  0.0001487756244077903  5.961041556248072e-7  9.255272505907766e-7
 2048  128  1  128  0  1.0859089885696814e-6  1.085908988611323e-6  1.0859089885993595e-6
 2048  128  1  128  1  3.232119967212712e-9  3.232120189218084e-9  3.232120190170034e-9
────────────────────────────────────────────────────────────────────────────────────────

As we see, the one-level $(j=1)$ additional correction bases from the eho-LOD method pushes the error obtained in the p-LOD/sp-LOD method down by an order of magnitude. The effects of the stabilization plus the additional corrections also compound, meaning that the stabilized bases along with the optimal correction-level $j$, yields the best solution.

Notes

1. wave_equation.jl writes a bunch of text files (something like ref_sol_n[FINE_SCALE]_T[TIME_STAMP].txt). In order to use the reference solution after generating once, use the option -r or --reference_sol in the command line.
2. If the output needs to be written to a file, use the -o or --output_file option while running the script from the command line.
3. The procedure to run the 2D examples located in the pLOD2d folder is identical to the 1D case. A few quick examples:

cd pLOD2d/
julia --project=. poisson.jl -n 128 -N 16 -p 1 -l 2 
julia --project=. wave_equation.jl -n 128 -N 4 -p 3 -l 4 -j 0
julia --project=. wave_equation.jl -n 128 -N 4 -p 3 -l 4 -j 1
julia --project=. wave_equation.jl -n 128 -N 4 -p 3 -l 4 -j 2

(Poisson Equation)
─────────────────────────────────────────────────────────────────────────────────
 1/h  1/H  p  l                pLOD          spLOD, HLM26           spLOD, DHM23 
─────────────────────────────────────────────────────────────────────────────────
 128   16  1  2  0.0186023684187968  0.001021492258520529  0.0026493359559217277
─────────────────────────────────────────────────────────────────────────────────

(Wave Equation)
─────────────────────────────────────────
 1/h  1/H  p  l  j          spLOD, HLM26 
─────────────────────────────────────────
 128    4  3  4  0  2.7283281596860636e-5
 128    4  3  4  1  4.935513468747671e-8
 128    4  3  4  2  3.3902371379255095e-8
──────────────────────────────────────────

References

1. Maier, R. 2021, SIAM Journal on Numerical Analysis 59(2), 1067-1089.
2. Dong, Z., Hauck, M., & Maier, R. (2023), SIAM Journal on Numerical Analysis, 61(4), 1918–1937.
3. Hauck M., Lozinski A., & Maier R. (2026) ESAIM Mathematical Modelling and Numerical Analysis, 60(1), 445–471.
4. Kalyanaraman, B., Krumbiegel, F., Maier, R., & Wang, S. (2025), arXiv [Math.NA].
5. Kalyanaraman, B., Krumbiegel, F., Maier, R., & Wang, S. (2026), arXiv [Math.NA].