Iterative parallel computing with PETSc

FreeFEM installation
1

Hi FreeFem++ developers!

Thank you very much for using FF.
I’d like to repeat the parallel computation using PETSc for implementing topology optimization.
So, I begin with the example elasticity-3d_PETSc.edp from FreeFem examples.

At first, I repeated the same analysis with two different meshes, run with the following code;

//Begin Program

load "medit"
load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

// Define and Load mesh
mesh3		Sh=readmesh3("./ShTrunced_25");
mesh3		ShPlt=Sh;
mesh3 		ShPltMoved;
plot(Sh,cmm="FirstMesh");

// Define Boundary Label
int[int] wall(3);     		wall    	= [1,7,13];			// 	fixed displacement
int[int] traction(1);     	traction    = [9];				//	traction

// 
int[int] n2o;
macro ShN2O()n2o// EOM
func Pk = [P2, P2, P2]; // finite element space
fespace Wh(Sh, Pk);
fespace WhPlt(ShPlt, Pk);

// Parameters
real Young = 1.0e8;
real poisson = 0.45;
real tmp = 1.0 + poisson;
real mu = Young  / (2.0 * tmp);
real lambda = Young * poisson / (tmp * (1.0 - 2.0 * poisson));

// Define Solid Mechanics
macro def(i)[i, i#B, i#C]// EOM     // vector field definition
macro init(i)[i, i, i]// EOM        // vector field initialization
macro gu 	[gu1,gu2,gu3] 				// Glabal	displacement vector
macro epsilon(u) 	[dx(u),dy(u#B),dz(u#C),(dz(u#B)+dy(u#C)), (dz(u)+dx(u#C)), (dy(u)+dx(u#B))] // strain tensor
macro D [
			[2.*mu+lambda,	lambda,			lambda,			0,				0,				0			], 
			[lambda,		2.*mu+lambda,	lambda,			0,				0,				0			], 
			[lambda,		lambda,			2.*mu+lambda,	0,				0,				0			], 
			[0,				0,				0,				2.*mu+lambda,	0,				0			], 
			[0,				0,				0,				0,				2.*mu+lambda,	0			], 
			[0,				0,				0,				0,				0,				2.*mu+lambda]
		] //elastic tensor
//macro g 	[0,0,-1000000.0] 				// traction vector
real gZ;	gZ = -1000000.0;
macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM

WhPlt [pltx, plty, pltz];
WhPlt [reducex, reducey, reducez];
WhPlt [gu1, gu2, gu3];
Wh def(u);
real[int] sol;
int[int] rest;

varf vPb(def(u), def(v)) = int3d(Sh)((D*epsilon(u))'*epsilon(v))
						+int2d(Sh,traction)(gZ*vC)
						+ on(wall, u = 0.0, uB = 0.0, uC = 0.0); // '
Mat A;


//---------------------------------1stAnalysis---------------------------------

buildDmesh(Sh);
createMat(Sh, A, Pk)

A = vPb(Wh, Wh);
real[int] rhs = vPb(0, Wh);
set(A, sparams = "-pc_type cholesky -ksp_view -ksp_max_it 100", bs = 3);
u[] = A^-1 * rhs;

if(!NoGraphicWindow) {
    ChangeNumbering(A, u[], sol);
    ChangeNumbering(A, u[], sol, inverse = true);
    rest = restrict(Wh, WhPlt, n2o);
    for[i, v : rest] pltx[][v] = u[][i];
    mpiReduce(pltx[], reducex[], processor(0, mpiCommWorld), mpiSUM);
}
gu1[]=reducex[];
gu2[]=reducey[];
gu3[]=reducez[];
ShPltMoved = movemesh3(ShPlt,transfo=[x+gu1,y+gu2,z+gu3]);
plot(ShPltMoved,cmm="Deformation1");

//-----------------------------------------------------------------------------

//Load Next Different Mesh
Sh=readmesh3("./ShTrunced_50");
ShPlt=Sh;
plot(Sh,cmm="SecondMesh");

//---------------------------------2ndAnalysis----------------------------------

buildDmesh(Sh);
createMat(Sh, A, Pk)

A = vPb(Wh, Wh);
rhs = vPb(0, Wh);
set(A, sparams = "-pc_type cholesky -ksp_view -ksp_max_it 100", bs = 3);
u[] = A^-1 * rhs;

if(!NoGraphicWindow) {
    ChangeNumbering(A, u[], sol);
    ChangeNumbering(A, u[], sol, inverse = true);
    rest = restrict(Wh, WhPlt, n2o);
    for[i, v : rest] pltx[][v] = u[][i];
    mpiReduce(pltx[], reducex[], processor(0, mpiCommWorld), mpiSUM);
}
gu1[]=reducex[];
gu2[]=reducey[];
gu3[]=reducez[];
ShPltMoved = movemesh3(ShPlt,transfo=[x+gu1,y+gu2,z+gu3]);
plot(ShPltMoved,cmm="Deformation2");

//End Program

When I do that, I get an error.

current line = 96 mpirank 3 / 4
current line = 96 mpirank 0 / 4
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** current line = 96 mpirank 2 / 4**
** current line = 96 mpirank 1 / 4**
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** err code 8 , mpirank 0**
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** err code 8 , mpirank 3**
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** err code 8 , mpirank 1**
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** err code 8 , mpirank 2**

I confirmed that multiple analysis runs can be performed on the same mesh.
I don’t know what procedures are required in case of a mesh change.
I’m very sorry if a similar question has already been asked and this is a rudimentary question.
I am sure you are very busy, but I would appreciate your guidance.

2

We cannot copy/paste your code. Make sure you use preformatted text.

3

Thank you very much for your confirmation and reply.
I apologize for any inconvenience caused. Thank you for pointing this out.

The text has been corrected to preformed text.
But here’s another thing I noticed.
To turn the above code, a mesh file is required but cannot be uploaded.
When I tried to upload it, I got the following response

Sorry, new users can not upload attachments.

Can you give me uploading privileges?

4

The error message is quite self-explanatory: Exec error : AssembleLinearForm size rhs and nb of DF of Vh. You need to resize your array rhs once you have loaded the new mesh.

5

Apologies if this is a silly question, but how do you resize it? Is there a line of code that you’d recommend adding?

6

Just .resize(), with the proper dimension.

7

Thanks, I’ll try this. I’ll probably put a full post up tomorrow.