// run with MPI: ff-mpirun -np 4 script.edp // NBPROC 4 cout.precision(17); if(mpirank==0){ cout.precision(17); macro grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient func Pk = P1; // finite element space mesh Th2 = square(40,40); // global mesh fespace Wh2(Th2, Pk); // local finite element space Wh2 du; varf varf1(du, v) = int2d(Th2)(grad(du)' * grad(v)+du*v) + on(1, du = 0.1 ); varf varf2(dum, v) = int2d(Th2)(grad(du)' * grad(v)+du*v) + on(1, dum = 0.1); matrix A2 = varf1(Wh2, Wh2,tgv = -1); real[int] rhs = varf1(0, Wh2, tgv = -10); du[] = A2^-1*rhs; real[int] out0 = A2*du[]; real[int] out1 = varf1(0,Wh2,tgv=-1); real[int] out2 = varf2(0,Wh2,tgv=-1); real res0 = out0'*out0; real res1 = out1'*out1; real res2 = out2'*out2; cout << " CASE 1: serial (w/o PETSc)" << endl; cout << "\tvalue 0 = " << res0 << endl; cout << "\tvalue 1 = " << res1 << endl; cout << "\tvalue 2 = " << res2 << endl; cout << " All values are equal. " << endl; } load "PETSc" // PETSc plugin macro dimension()2// EOM // 2D or 3D include "macro_ddm.idp" // additional DDM functions macro grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient func Pk = P1; // finite element space mesh Th = square(getARGV("-global", 40), getARGV("-global", 40)); // global mesh Mat A; buildMat(Th, getARGV("-split", 1), A, Pk, mpiCommWorld) fespace Wh(Th, Pk); // local finite element space Wh du; varf varf1(du, v) = int2d(Th)(grad(du)' * grad(v)+du*v) + on(1, du = 0.1 ); varf varf2(dum, v) = int2d(Th)(grad(du)' * grad(v)+du*v) + on(1, dum = 0.1); A = varf1(Wh, Wh,tgv = -1); real[int] rhs = varf1(0, Wh, tgv = -10); du[] = A^-1*rhs; real[int] out0 = A*du[]; real[int] out1 = varf1(0,Wh,tgv=-1); real[int] out2 = varf2(0,Wh,tgv=-1); real res0 = A(out0,out0); real res1 = A(out1,out1); real res2 = A(out2,out2); if(mpirank==0){ cout << " CASE 2: with PETSc" << endl; cout << "\tvalue 0 = " << res0 << endl; cout << "\tvalue 1 = " << res1 << endl; cout << "\tvalue 2 = " << res2 << endl; cout << " I understand why value 0 and value 1 may not be equal... " << endl; cout << " But why does value 1 not equal value 2 in case 2 (even with 1 proc)? " << endl; }