Integral over nested meshes

General Discussion
1

Dear all,

I am trying to evaluate an integral for the solution obtained in the nested meshes.
More concretely, I’d like to evaluate the integral of the solution in the following example:

According to the previous entry (Issue with transfer between two meshes),
Maybe, I should add the n2o operator for the mesh “ThFine”, as in the example “maxwell-3d-PETSc.edp”.
However, I don’t know how to do this with “buildMatEdgeRecursive.”

Would someone show me how to do this?
Or, is there any other method?

Thank you in advance.

2

Instead of computing an integral, you want to compute a weighted integral, using the partition of unity.
E.g., instead of

real val = int3d(Th)(real(u));
mpiAllReduce(val, ..., mpiSUM);

Do

real[int] D = MG[0].D; // partition of unity on the fine mesh
u .*= D; // entry-wise scaling by the partition of unity
real val = int3d(Th)(u);
mpiAllReduce(val, ..., mpiSUM);
3

Dear prj,
Thank you so much for your response.

To confirm the integral value, I evaluated the volume of a unit cube based on your method.

//volume of a unit cube
fespace VhP1(ThFine, P1);
VhP1 volint = 1.0;
real valwrong = int3d(ThFine)(volint);
real valwrongG;
mpiAllReduce(valwrong, valwrongG, mpiCommWorld, mpiSUM);
real[int] D = MG[0].D; // partition of unity on the fine mesh, OK
volint[].*=D; // entry-wise scaling by the partition of unity
real val = int3d(ThFine)(volint);
real valG;
mpiAllReduce(val, valG, mpiCommWorld, mpiSUM);
if(mpirank==0){
  cout << "valwrongG = " << valwrongG << endl;
  cout << "valG = " << valG << endl;
}

If I run the code with np=1, I obtained correct results with valwrongG=1 and valG=1.
However, when I run the code with np=4, they are evaluated as valwrongG=1.709913 and valG = 8.724715e-01. That is, both of them are apart from 1, and even valG seems to contain an error.

Am I missing something?

4

OK, there are multiple reasons why it’s not working (one is that you are mixing the partition of unity which is of Nédélec type with a P_1 function). Still, there are some issues in macro_ddm.idp, or missing pieces, to make this work. I’ll try to think of something and get back to you ASAP. Sorry about that.

5

I have a solution, far from ideal, but it gets the job done. Please fetch Adds buildMatEdgeRecursiveWithPartitioning. · FreeFem/FreeFem-sources@57d5ad4 · GitHub.
Then, the following code works and always returns a volume of 1.

diff --git a/examples/hpddm/maxwell-mg-3d-PETSc-complex.edp b/examples/hpddm/maxwell-mg-3d-PETSc-complex.edp
index a70f15b8..c1a872ba 100644
--- a/examples/hpddm/maxwell-mg-3d-PETSc-complex.edp
+++ b/examples/hpddm/maxwell-mg-3d-PETSc-complex.edp
@@ -28,3 +28,25 @@ func Pk = Edge03d;
 func PkPart = Edge03ds0;
-buildMatEdgeRecursive(Th, s, level, P, MG, Pk, mpiCommWorld, PkPart, defPart, initPart);
+fespace PhFine(Th[0], P0);
+PhFine partFine;
+{
+    fespace PhCoarse(Th[level - 1], P0);
+    PhCoarse part;
+    partitionerSeq(part[], Th[level - 1], mpisize);
+    partitionerPar(part[], Th[level - 1], mpiCommSelf, mpisize);
+    buildMatEdgeRecursiveWithPartitioning(Th, part[], s, level, P, MG, Pk, mpiCommWorld, PkPart, defPart, initPart);
+    part = part;
+    partFine = part;
+}
+
+{
+    // volume of a unit cube
+    mesh3 ThNo = trunc(Th[0], abs(partFine-mpirank) < 1e-2);
+    fespace VhP1(ThNo, P1);
+    VhP1 volint = 1.0;
+    real val = int3d(ThNo)(volint);
+    real valG;
+    mpiAllReduce(val, valG, mpiCommWorld, mpiSUM);
+    if(mpirank==0)
+        cout << "valG = " << valG << endl;
+}

The idea is to let you store the partitioning, so that you can get a nonoverlapping decomposition on the fine level, so that you can integrate properly any variable without duplicated unknowns. Let me know if something is not clear to you.

6

Thank you so much! It worked very well in my environment.

7

Dear prj,

Thank you so much for this solution.
I have an additional question.
I would like to solve Maxwell’s equation with periodic boundary conditions based on the example in this post.
However, the number of nodes on the surfaces of a cube does not match if I use your macro “buildMatEdgeRecursiveWithPartitioning” because of the partitions.
I think a macro like
macro ThPeriodicity()labPeriodic//
will solve this problem as in the following example.

However, I don’t know how to define such a macro because we define the mesh with an array like
Th[level - 1] = cube(n, n, n, [x, y, z], label = LL);

It would be very helpful if you show me the solution when you are available.
Thank you in advance.

8

I believe they should match. Could you please provide a full running example that highlight the exact issue?

9

Dear prj,

Thank you for your response.
I modified the example a little bit as follows:

//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4

load "PETSc-complex"
load "Element_Mixte3d"
macro dimension()3// EOM
include "macro_ddm.idp"
include "cube.idp"

macro Curl(ux, uy, uz)[dy(uz)-dz(uy), dz(ux)-dx(uz), dx(uy)-dy(ux)]// EOM
macro CrossN(ux, uy, uz)[uy*N.z-uz*N.y, uz*N.x-ux*N.z, ux*N.y-uy*N.x]// EOM
macro def(i)[i, i#y, i#z]// EOM
macro init(i)[i, i, i]// EOM
macro defPart(i)i// EOM
macro initPart(i)i// EOM

int level = 2;
int s = 3;

Mat<complex>[int] MG(level);
mesh3[int] Th(level);
matrix[int] P(level - 1);
//int Robin = 2;
//int[int] LL = [Robin, Robin, Robin, Robin, Robin, Robin];
real n = getARGV("-n", 21) / s;
Th[level - 1] = cube(n, n, n, [x, y, z]);
func perio = [[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]];

//func Pk = Edge03d;
//func PkPart = Edge03ds0;
macro Pk() Edge03d, periodic = perio//EOM
macro PkPart() Edge03ds0, periodic = perio//EOM

//buildMatEdgeRecursive(Th, s, level, P, MG, Pk, mpiCommWorld, PkPart, defPart, initPart);
fespace PhFine(Th[0], P0);
PhFine partFine;
{
    fespace PhCoarse(Th[level - 1], P0);
    PhCoarse part;
    partitionerSeq(part[], Th[level - 1], mpisize);
    partitionerPar(part[], Th[level - 1], mpiCommSelf, mpisize);
    buildMatEdgeRecursiveWithPartitioning(Th, part[], s, level, P, MG, Pk, mpiCommWorld, PkPart, defPart, initPart);
    part = part;
    partFine = part;
}

func k = 6 * pi;
func f = exp(-60 * ((x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2));

complex[int] rhs;
matrix<complex>[int] Opt(level);
for(int i = 0; i < level; ++i) {
    fespace Vh(Th[i], Pk);
    varf vPb([Ex,Ey,Ez], [vx,vy,vz]) =
        int3d(Th[i])(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
      + int3d(Th[i])(-k^2*[vx,vy,vz]'*[Ex,Ey,Ez])
      //+ int2d(Th[i],Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
      + int2d(Th[i],-111111)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez));
    Opt[i] = vPb(Vh, Vh, sym = 1);
    MG[i] = Opt[i];
    if(i == 0) {
        varf vRHS([Ex,Ey,Ez],[vx,vy,vz]) = -int3d(Th[i])([vx,vy,vz]'*[0,0,f]);
        rhs.resize(Vh.ndof);
        rhs = vRHS(0, Vh);
    }
}
set(MG, P, sparams = "-pc_type mg -ksp_monitor -ksp_view_final_residual -ksp_type fgmres -ksp_gmres_restart 200 -ksp_max_it 200");
set(MG, 0, sparams = "-mg_coarse_ksp_type gmres -mg_coarse_ksp_rtol 1e-1 -mg_coarse_ksp_pc_side right " + " -mg_coarse_ksp_max_it 50 -mg_coarse_ksp_converged_reason -mg_coarse_pc_type asm -mg_coarse_sub_pc_factor_mat_solver_type mumps -mg_coarse_sub_pc_type cholesky -mg_coarse_pc_asm_type restrict", O = Opt[level - 1]);
set(MG, level - 1, sparams = "-mg_levels_ksp_type richardson -mg_levels_ksp_pc_side left -mg_levels_pc_type asm -mg_levels_sub_pc_type cholesky -mg_levels_sub_pc_factor_mat_solver_type mumps -mg_levels_pc_asm_type restrict", O = Opt[0]);
fespace Vh(Th[0], Pk);
Vh<complex> def(u);
u[] = MG[0]^-1 * rhs;
mesh3 ThFine = Th[0];
fespace VhFine(ThFine, P1);
VhFine plt = sqrt(real(u)^2 + real(uy)^2 + real(uz)^2);
plotD(ThFine, plt, cmm = "Solution");
int[int] fforder = [1, 1];
savevtk("maxwell-mg-3d.vtu", ThFine, [real(u), real(uy), real(uz)], [imag(u), imag(uy), imag(uz)], order = fforder, bin = 1);

It runs with np=1. However, an error about the number of nodes on the surfaces appears when n>=2.

10

This is not expected to work, you need to make sure your partitioning will preserve periodicity, see, e.g., FreeFem-sources/diffusion-periodic-balanced-2d-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub.

11

Thank you for your reply.
In my understanding, I need to combine the codes in lines 24–105 in “diffusion-periodic-balanced-2d-PETSc.edp” with your macro “buildMatEdgeRecursiveWithPartitioning.” Is this right?
I will modify the code based on your advice.

12

There are many ways to achieve what you need. Another way is to decompose the initial coarse cube in two, a part touching the periodic boundary, and the rest. Then call partitionerSeq() on the rest. The partitioning made up of “part touching the periodic boundary” + “the rest + 1” will preserve periodicity. Process #0 will be in charge of the part touching the periodic boundary, and other processes will handle the rest of the domain.

13

I understood the concept well. I would like to implement the periodic boundary conditions using one of the methods you suggested. Thank you so much!

14

Let me know if you’re stuck, I’d be happy to help.