Can I run a discontinous galerkin code for periodic boundary conditions in free fem ++? Which version is good for this?
Yes but it is written in the formulation due to discontinuity.
Can you please explain how to write the boundary terms in the code?
Which DG formulation do you use ?
This a math problem not a freefem problem , but the idea is simple,
you periodic Boundary condition because like internal border so with the same scheme as for internal edge, but the jump and mean must be computed by hand fo exemple
jump(u)*jump(v) on border 2 and 4 become
Gamma_2 is the translation Tp of Gamma_4 : (x,y) → (x-1,y)
int1d(Th,2) ( (u-u(Tp)) ( v - v(Tp)) now the problem is to compute this integral
remenber you can move the quadrature point so you can slip the formula in 4 term
( no tested so may by error!)
int1d(Th,2) ( (u) ( v ))
+ int1d(Th,2,matu=Tp) ( - (u(Tp) ( v ))
+ int1d(Th,2,mapt=Tp) ( - (u) ( v(Tp) )
+ int1d(Th,4) ( (u v )
As you can see it is hard to do !
(post deleted by author)
Dear Professor
I try to write the jump term and flux term for the periodic boundary condition by using the mapu and mapt.
But the output is
Exec error : no map in the case (2)??
– number :1
It seems that correctly completing the periodic boundary conditions with DG in Freefem is not an easy thing.
Could you please explain it?
// solve poisson equation
macro flux(u, v) (mean(dx(u) * N.x + dy(u) * N.y) * jump(v) + mean(dx(v) * N.x + dy(v) * N.y) * jump(u)) // EOM
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), flags=1);
fespace Vh(Th, P1dc, periodic = [[2, y], [4, y]]);
Vh u, v;
solve Poisson(u, v) = int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v))
- intalledges(Th)((1 - nTonEdge) * flux(u,v) / nTonEdge)
+ intalledges(Th)((1 - nTonEdge) * (penalty / lenEdge * jump(u) * jump(v)) / nTonEdge)
// flux term
// jump term
+ int1d(Th, 2)((penalty / lenEdge * u * v) / 2)
- int1d(Th, 2, mapu=[x - 1.0, y])((penalty / lenEdge * u * v) / 2)
- int1d(Th, 2, mapt=[x - 1.0, y])((penalty / lenEdge * u * v) / 2)
+ int1d(Th, 4)((penalty / lenEdge * u * v) / 2)
- 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++){
int Nx = 2 ^ (i + 2);
cout << "1/" << Nx << ":\t";
cout.scientific << errL2(i) << "\t";
cout.fixed << orderL2(i) << "\t";
cout.scientific << errInf(i) << "\t";
cout.fixed << orderInf(i) << "\t";
cout.scientific << errH1(i) << "\t";
cout.fixed << orderH1(i) << "\t" << endl;
}
A bug fix has been applied on version 4.15
“fix problem in integration of moving test or unknown function in 2d mesh:”
Do you use version 4.15, or older version?
Thanks for your reply!!
I just deleted the old version of Freefem, downloaded the version 4.15 and ran this code.
The problem still occurs
– FreeFem++ v4.15 (Mer 11 déc 2024 13:31:45 CET - git v4.15)
file : poisson_dg_SIPG_test.edp
Load: lg_fem lg_mesh lg_mesh3 eigenvalue
1 : // solve poisson equation
2 : macro flux(u, v) (mean(dx(u) * N.x + dy(u) * N.y) * jump(v) + mean(dx(v) * N.x + dy(v) * N.y) * jump(u)) ) // EOM
3 : int MaxIter = 6;
4 : // mesh size [4 * 4], [8 * 8], … ,[2 ^ (MaxIter + 1), 2 ^ (MaxIter + 1)]
5 : real[int] errL2(MaxIter), errInf(MaxIter), errH1(MaxIter);
6 : real[int] orderL2(MaxIter), orderInf(MaxIter), orderH1(MaxIter);
7 : func rhs = 8 * pi ^ 2 * sin(2 * pi * x) * cos(2 * pi * y);
8 : func uExact = sin(2 * pi * x) * cos(2 * pi * y);
9 : func udxExact = 2 * pi * cos(2 * pi * x) * cos(2 * pi * y);
10 : func udyExact = -2 * pi * sin(2 * pi * x) * sin(2 * pi * y);
11 : func uBcd = 0;
12 : real cpu = clock();
13 : real penalty = 20.0;
14 : // mesh iteration
15 : for (int i = 0; i < MaxIter; i++){
16 : mesh Th = square(2 ^ (i + 2), 2 ^ (i + 2), flags=1);
17 : fespace Vh(Th, P1dc, periodic = [[2, y], [4, y]]);
18 : Vh u, v;
19 : solve Poisson(u, v) = int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v))
20 : - intalledges(Th)((1 - nTonEdge) * flux(u,v) (mean(dx(u) * N.x + dy(u) * N.y) * jump(v) + mean(dx(v) * N.x + dy(v) * N.y) * jump(u)) / nTonEdge)
21 : + intalledges(Th)((1 - nTonEdge) * (penalty / lenEdge * jump(u) * jump(v)) / nTonEdge)
22 : // flux term
23 :
24 : // jump term
25 : + int1d(Th, 2)((penalty / lenEdge * u * v) / 2)
26 : - int1d(Th, 2, mapu=[x - 1.0, y])((penalty / lenEdge * u * v) / 2)
27 : - int1d(Th, 2, mapt=[x - 1.0, y])((penalty / lenEdge * u * v) / 2)
28 : + int1d(Th, 4)((penalty / lenEdge * u * v) / 2)
29 : - int2d(Th)(rhs * v);
30 : // error
31 : errL2[i] = sqrt(int2d(Th)((u - uExact) ^ 2));
32 : Vh absError = abs(u - uExact);
33 : errInf[i] = absError.max;
34 : Vh udx = dx(u), udy = dy(u);
35 : errH1[i] = sqrt(int2d(Th)((udx - udxExact) ^ 2) + int2d(Th)((udy - udyExact) ^ 2));
36 : }
37 : for (int i = 1; i < MaxIter; i++){
38 : orderL2[i] = log(errL2[i - 1] / errL2[i]) / log(2.);
39 : orderInf[i] = log(errInf[i - 1] / errInf[i]) / log(2.);
40 : orderH1[i] = log(errH1[i - 1] / errH1[i]) / log(2.);
41 : }
42 : cout << “================OUTPUT===============” << endl;
43 : cout << "CPU time = " << clock() - cpu << endl;
44 : cout << "errors and convergence rates = " << endl;
45 : cout << “h L2 error order Inf error order H1_semi e
… : rror order” << endl;
46 : for (int i = 0; i < MaxIter; i++){
47 : int Nx = 2 ^ (i + 2);
48 : cout << “1/” << Nx << ": ";
49 : cout.scientific << errL2(i) << " ";
50 : cout.fixed << orderL2(i) << " ";
51 : cout.scientific << errInf(i) << " ";
52 : cout.fixed << orderInf(i) << " ";
53 : cout.scientific << errH1(i) << " ";
54 : cout.fixed << orderH1(i) << " " << endl;
55 : } sizestack + 1024 =9336 ( 8312 )– Square mesh : nb vertices =25 , nb triangles = 32 , nb boundary edges 16 rmdup= 0
current line = 29
Exec error : no map in the case (2)??
– number :1
Exec error : no map in the case (2)??
– number :1
err code 8 , mpirank 0
Then I don’t know why this error occurs…