Bug with Petsc in building a local matrix with edge elements

General Discussion
1

Dear all,

I encounter a bug when assembling a matrix with Nedelec elements. The problem does not show up with 4 procs but with 8 and more.

Thank you for your help

//=================================================
load “msh3”
load “medit”
load “PETSc-complex”
//load “PETSc”
load “mmg”
load “tetgen”
include “MeshSurface.idp”
macro dimension()3// EOM // 2D or 3D
include “macro_ddm.idp” // additional DDM functions

//==================================================
// parameters (physical, numerics)
//==================================================
int R = 5;
real kappa = 4.;
real nodePerXi = 1;
int nbsegul = kappa*nodePerXi;
real hsize= 1./nbsegul;

// applied field
real Bax = 0.;
real Bay = 0.;
real Baz = 0.9;

real omega = 1.;
tgv = 1e30;
real dt = 0.1;
real idt=1./dt;

//==================================================
// files for savings
//==================================================

ofstream fout(“./a_echo”);
fout.precision(12);

//==================================================
// macros
//==================================================
macro div(A)(dx(A) + dy(A#y) + dz(A#z))//EOM
macro Grad(u)([dx(u), dy(u), dz(u)])//

//==================================================
//finite element space
//==================================================

func PkMagEdge = Edge03d;
func PkMagRT = RT03d;

//===============================================================
// global mesh
//===============================================================

//---------- sphere de référence
meshS Thsph0 = Sphere(R,hsize,1,1);
real[int] domain = [0.,0.,0.,1,hsize^3/6];
mesh3 Th=tetg(Thsph0,switch=“paAAQYY”,nbofregions=1,regionlist=domain);
if(mpirank==0) medit(“sphere”, Th);

// if(mpirank==0)
// {
// cout << " ========== sphere.hmin = " << Thsph.hmin << endl;
// cout << " ========== sphere.hmax = " << Thsph.hmax << endl;
// cout << " ========== sphere.nt = " << Thsph.nt << endl;
// fout << " ========== sphere.hmin = " << Thsph.hmin << endl;
// fout << " ========== sphere.hmax = " << Thsph.hmax << endl;
// fout << " ========== sphere.nt = " << Thsph.nt << endl;
// }

//===============================================================
//local FE spaces
//===============================================================

buildDmesh(Th)

//local FE spaces
fespace VhMagEdge(Th, PkMagEdge);
fespace VhMagRT(Th, PkMagRT);

if(mpirank==0)
{
cout << " ********** local mesh on proc = " << mpirank << endl;
cout << " ========== ntetra = " << Th.nt << endl;
cout << " ========== dof magField = " << VhMagEdge.ndof << endl;
cout << " ========== dof vectPot = " << VhMagRT.ndof << endl;

fout << " ********** local mesh on proc = " << mpirank << endl;
fout << " ========== ntetra = " << Th.nt << endl;
fout << " ========== dof magField = " << VhMagEdge.ndof << endl;
fout << " ========== dof vectPot = " << VhMagRT.ndof << endl;

}

//===============================================================
//local matrices
//local weak forms
//===============================================================

//local matrices
Mat Asys11, Asys22;

{
macro def(i)[i, i#y, i#z]//EOM // vector field definition
macro init(i)[i, i, i]//EOM // vector field initialization
createMat(Th, Asys11, PkMagEdge);

// macro necessaire pour la construction de Asys22 (bug PETSc)
macro ThPostProcessD(D) {
VhMagRT def(u), def(v);
varf onG(def(u), def(v)) = on(-111111, u = 10 + x - z, uy = 100 + y - x, uz = 1000 + z - y);
v = D;
u = onG(0, VhMagRT);
for [j, dj : u] dj = abs(dj) > 1e-2 ? 0.0 : 1.0;
def(u) = [u, uy, uz];
D = u;
} //EOM
createMat(Th, Asys22, PkMagRT);
}

// local weak forms

varf vSys11([sigma, sigmay, sigmaz],[ki, kiy, kiz]) =
int3d(Th, qfV = qfV2)([sigma, sigmay, sigmaz]'*[ki, kiy, kiz])

  • on(1, sigma = Bax, sigmay = Bay, sigmaz = Baz);

varf vSys22([A, Ay, Az],[v, vy, vz]) =
int3d(Th, qfV = qfV2)(idt*[A, Ay, Az]'*[v, vy, vz])

  • int3d(Th, qfV = qfV2)(omega*div(A)*div(v))
  • on(1, A = 0., Ay = 0., Az = 0.); // [A, Ay, Az] = 0 ==> A.n = 0 here

// local matrices
Asys11 = vSys11(VhMagEdge, VhMagEdge); // bug ici avec 16 proc mais 4 proc ok
//Asys22 = vSys22(VhMagRT, VhMagRT);

2

Can’t copy/paste your code => can’t help you.

3

Here is the code:

(attachments)

Zcode.edp (4.04 KB)

4

I see two wrong things:

  1. you should broadcast your mesh from processor 0 after the mesh adaptation step because there are no guarantees that the mesh will be the same among all processes
  2. for edge elements, you should not define init and part but instead PkPart, see FreeFem-sources/examples/hpddm/maxwell-3d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub

With both fixes, the code runs to completion:

$ ff-mpirun -n 16 Zcode.edp -v 0 -version
[...]
********** local mesh on proc = 0
 ========== ntetra = 29545
 ========== dof magField = 37680
 ========== dof vectPot = 61270
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There are 2 unused database options. They are:
Option left: name:-nw value: Zcode.edp source: command line
Option left: name:-v value: 0 source: command line
5

Thank you very much for your help