DG method with periodic boundary condition

kail · May 27, 2025, 12:11pm

I don’t know how to deal with the periodic bcd in DG method.
I have successfully finished the Dirichlet or Nuemann bcd for Poisson equation with DG method.
I set the left and right bcd is periodic, and the main difficult is how to deal jump and penalty term on left or right bcd.
Looking forward to one’s reply and help!

// solve poisson equation
int MaxIter = 6;
// mesh size [4 * 4], [8 * 8], ... ,[2 ^ (MaxIter + 1), 2 ^ (MaxIter + 1)]
real[int] errL2(MaxIter), errInf(MaxIter), errH1(MaxIter);
real[int] orderL2(MaxIter), orderInf(MaxIter), orderH1(MaxIter);
func rhs = 8 * pi ^ 2 * sin(2 * pi * x) * cos(2 * pi * y);
func uExact = sin(2 * pi * x) * cos(2 * pi * y);
func udxExact = 2 * pi * cos(2 * pi * x) * cos(2 * pi * y);
func udyExact = -2 * pi * sin(2 * pi * x) * sin(2 * pi * y);
func uBcd = 0;
real cpu = clock();
real penalty = 20.0; 
// mesh iteration
for (int i = 0; i < MaxIter; i++){
    mesh Th = square(2 ^ (i + 2), 2 ^ (i + 2));
    fespace Vh(Th, P1dc, periodic = [[2, y], [4, y]]);
    Vh u, v;
    solve Poisson(u, v, solver=UMFPACK) = int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v))
	- intalledges(Th)((1 - nTonEdge) * (mean(dx(u)) * jump(v) * N.x + mean(dy(u)) * jump(v) * N.y
        + mean(dx(v)) * jump(u) * N.x + mean(dy(v)) * jump(u) * N.y) / 2)
	+ intalledges(Th)((1 - nTonEdge) * (penalty / lenEdge * jump(u) * jump(v)) / 2)
    ////////////////////////////////////////////////////////////////////////////////////////////////////
	- intalledges(Th, 4)((mean(dx(u)) * jump(v) * N.x + mean(dy(u)) * jump(v) * N.y
        + mean(dx(v)) * jump(u) * N.x + mean(dy(v)) * jump(u) * N.y))
	+ intalledges(Th, 4)((penalty / lenEdge * jump(u) * jump(v)))
    ////////////////////////////////////////////////////////////////////////////////////////////////////
    - int2d(Th)(rhs * v);
    // error
    errL2[i] = sqrt(int2d(Th)((u - uExact) ^ 2));
    Vh absError = abs(u - uExact);
    errInf[i] = absError[].max;
    Vh udx = dx(u), udy = dy(u);
    errH1[i] = sqrt(int2d(Th)((udx - udxExact) ^ 2) + int2d(Th)((udy - udyExact) ^ 2));
}
for (int i = 1; i < MaxIter; i++){
    orderL2[i] = log(errL2[i - 1] / errL2[i]) / log(2.);
    orderInf[i] = log(errInf[i - 1] / errInf[i]) / log(2.);
    orderH1[i] = log(errH1[i - 1] / errH1[i]) / log(2.);
}
cout << "================OUTPUT===============" << endl;
cout << "CPU time = " << clock() - cpu << endl;
cout << "errors and convergence rates = " << endl;
cout << "h\t L2 error\t order\t\t Inf error\t order\t\t H1_semi error\t order" << endl;
for (int i = 0; i < MaxIter; i++){
    cout << "1/" << 2 ^ (i + 2) << ":\t" 
    << errL2[i] << "\t" << orderL2[i] << "\t\t" 
    << errInf[i] << "\t" << orderInf[i] << "\t\t" 
    << errH1[i] << "\t" << orderH1[i] << "\t" << endl;
}

kail · May 29, 2025, 1:57am

Dear all
After a day of trying, I still couldn’t write the correct code. I finally decided to give up this problem, hoping that it could be solved one day

// solve poisson equation
int MaxIter = 1;
// mesh size [4 * 4], [8 * 8], … ,[2 ^ (MaxIter + 1), 2 ^ (MaxIter + 1)]
real[int] errL2(MaxIter), errInf(MaxIter), errH1(MaxIter);
real[int] orderL2(MaxIter), orderInf(MaxIter), orderH1(MaxIter);
func rhs = 8 * pi ^ 2 * sin(2 * pi * x) * cos(2 * pi * y);
func uExact = sin(2 * pi * x) * cos(2 * pi * y);
func udxExact = 2 * pi * cos(2 * pi * x) * cos(2 * pi * y);
func udyExact = -2 * pi * sin(2 * pi * x) * sin(2 * pi * y);
func uBcd = 0;
// func real[int] Tp(real[int] X) {
// real[int] Y(2);
// Y[0] = X[0] - 1.0;
// Y[1] = X[1];
// return Y;
// }
real cpu = clock();
real penalty = 20.0;
// mesh iteration
for (int i = 0; i < MaxIter; i++){
int n = 2^(i+2);
mesh Th = square(2 ^ (i + 2), 2 ^ (i + 2), flags=1);
fespace Vh(Th, P1dc, periodic = [[2, y], [4, y]]);
Vh u, v;
problem Poisson(u, v, solver=UMFPACK) = int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v))
- intalledges(Th)((1 - nTonEdge) * (mean(dx(u)) * jump(v) * N.x + mean(dy(u)) * jump(v) * N.y
+ mean(dx(v)) * jump(u) * N.x + mean(dy(v)) * jump(u) * N.y) / 2)
+ intalledges(Th)((1 - nTonEdge) * (penalty / lenEdge * jump(u) * jump(v)) / 2)
//////////////////////////////////////////////////
// - int1d(Th, 2)((mean(dx(u)) * jump(v) * N.x + mean(dy(u)) * jump(v) * N.y
// + mean(dx(v)) * jump(u) * N.x + mean(dy(v)) * jump(u) * N.y))
// + int1d(Th, 2)((penalty / lenEdge * jump(u) * jump(v)))
+ int1d(Th, 2)(penalty / lenEdge * u * v)
- int1d(Th, 2, mapu=[x-1, y])(penalty / lenEdge * u * v)
- int1d(Th, 2, mapt=[x-1, y])(penalty / lenEdge * u * v)
+ int1d(Th, 4)(penalty / lenEdge * u * v)
//////////////////////////////////////////////////
- int2d(Th)(rhs * v);
Poisson;
// error
errL2[i] = sqrt(int2d(Th)((u - uExact) ^ 2));
Vh absError = abs(u - uExact);
errInf[i] = absError[].max;
Vh udx = dx(u), udy = dy(u);
errH1[i] = sqrt(int2d(Th)((udx - udxExact) ^ 2) + int2d(Th)((udy - udyExact) ^ 2));
}
for (int i = 1; i < MaxIter; i++){
orderL2[i] = log(errL2[i - 1] / errL2[i]) / log(2.);
orderInf[i] = log(errInf[i - 1] / errInf[i]) / log(2.);
orderH1[i] = log(errH1[i - 1] / errH1[i]) / log(2.);
}
cout << “================OUTPUT===============” << endl;
cout << "CPU time = " << clock() - cpu << endl;
cout << "errors and convergence rates = " << endl;
cout << “h\t L2 error\t order\t\t Inf error\t order\t\t H1_semi error\t order” << endl;
for (int i = 0; i < MaxIter; i++){
cout << “1/” << 2 ^ (i + 2) << “:\t”
<< errL2[i] << “\t” << orderL2[i] << “\t\t”
<< errInf[i] << “\t” << orderInf[i] << “\t\t”
<< errH1[i] << “\t” << orderH1[i] << “\t” << endl;
}

Monirul25 · May 29, 2025, 5:40am

Could you please upload your scheme here, along with all the details you intend to implement in the FreeFEM++ code?.

kail · May 29, 2025, 6:28am

Thanks for your reply!!

To present the problem more clearly, I uploaded the code of this equation under the Dirichlet and Nuemann boundary conditions.

I think the key challenge lies in integrating along the periodic boundary by treating it as an interior edge.
poisson_dg_SIPG_dirichlet.edp (2.2 KB)
poisson_dg_SIPG_nuemann.edp (2.1 KB)

Monirul25 · May 29, 2025, 6:37am

Could you please tell me what periodic boundary conditions you want to implement in the code?.

kail · May 29, 2025, 6:52am

Thank you so much for your response!
I want to implement the periodic boundary conditions on the left (4) and right (2) boundaries, i.e., u(0,y)=u(1,y).

fespace Vh(Th, P1dc, periodic = [[2, y], [4, y]]);

The upper and lower boundaries I use 0 nuemann.

Monirul25 · May 29, 2025, 6:54am

Poisson_Periodic_DG_MI.edp (3.5 KB)
(Can you please arrange the output properly).

I hope results will comes goods. There are some silly errors in codes. I will check it. Can you send me the code by arranging the output properly?.

(Need to write infimum error and H1_semi error properly format).

kail · May 29, 2025, 7:23am

I have modified the layout of the output
Poisson_Periodic_DG_MI.edp (3.1 KB)

Monirul25 · May 29, 2025, 7:34am

Poisson_Periodic_DG_MI.edp (4.6 KB)
(Increase your penality parameter).

Monirul25 · May 29, 2025, 7:39am

Poisson_Periodic_DG_MI(Latest).edp (3.1 KB)
(Result getting improve if you increase penality parameter)

kail · May 29, 2025, 7:46am

Thanks for your reply!
I observed that under the low penalty parameter, when I set the number of grids to 512 or bigger, the error was not consistent with the theoretical result.
Is this because the jump term was not considered in int1d(Th, 2) and int1d(Th, 4)?
Should this part be taken into account?

Monirul25 · May 29, 2025, 7:50am

Your intalleges terms are not properly wriiten. I will share a refine version of it to you today.
Regarding low penality, penality parameter depends on methods.

kail · May 29, 2025, 7:53am

OK! I am extremely grateful for your sharing and correction!

kail · May 30, 2025, 2:11am

Dear Monirul,

I’m very sorry to bother you again, but I’m still having trouble identifying what is wrong with the intalledges terms in my DG implementation.

I would greatly appreciate it if you could share your refined version of the intalledges terms.

Thank you so much for your time and help.

Monirul25 · May 30, 2025, 2:31am

Your code is correct but not written in a proper way. If you use definition of macro for the Neumann derivative term inside intalledges term then it will look goods. Okay. I will share you today.

kail · May 30, 2025, 2:45am

Sure, really appreciate it!!!

Monirul25 · May 30, 2025, 3:41am

Your formulation are different brother. Macro definition is not working in your formulations.
I have seen in intalledges term, you are dividing by 2. But it should be divide by nTonEdge although your code gives good results.

kail · May 30, 2025, 4:33am

nTonEdge = 2 on the inner edge, so I directly divide by 2.
Anyway, thanks a lot for your reply!

Topic		Replies	Views
dG code for periodic boundary condition General Discussion	3	77	April 3, 2025
Mixed Boundary Condition and Periodic System General Discussion	0	466	March 17, 2021
"Jump" boundary condition, or a Periodic Boundary Condition with offset General Discussion	1	568	April 15, 2022
Periodic boundary about taylor problem General Discussion	1	407	January 6, 2021
Periodic BC with DG spaces General Discussion	1	109	April 16, 2024

DG method with periodic boundary condition

Related topics