Include two different PDE problems in PETSc

General Discussion
1

Hi FreeFem++ developers!

Thank you very much for makingFF.

I’d like to perform the 3D topology optimization, where the 3D elasticity problem and reaction-diffusion equation (RDE) are included, using PETSc.

I am now attempting to add a parallelized RDE to my already parallelized 3D elasticity problem code as follows.

verbosity = 0;

load "medit" 	// for data save
load "msh3"
load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"  

real E;			E 			= 1.0;						// Young's modulus	
real nu; 		nu 			= 0.3;						// Poisson's ratio
real lambda;	lambda 		= nu*1./(1+nu)/(1-2.*nu);	// Lame coefficient
real mu; 		mu 			= 1./2./(1+nu);				// Lame coefficient
real matd;		matd 		= 1e-3;
real matw;		matw		= 0.5;

int[int] 	labs(6);		// labels on boundaries of squre: (bottom, right, top, left)
mesh3 		Shc1;				// Mesh of the squre 1
labs		= [1,0,3,0,0,6];
Shc1 		= cube(5,10,1,[-0.25+(0.25+0.25)*x, 0.0+(1.0-0.0)*y, -0.25+(-0.15+0.25)*z],label=labs,flags=1,region=1);
mesh3 		Shc2;				// Mesh of the squre  2
labs		= [7,8,9,100,11,12];	
Shc2 		= cube(5,10,4,[-0.25+(0.25+0.25)*x, 0.0+(1.0-0.0)*y, -0.15+(0.25+0.15)*z],label=labs,flags=1,region=1);
mesh3 		Sh=Shc1+Shc2;//+Shc3;
mesh3		ShPlt=Sh;
mesh3 		ShPltMoved;
mesh3 		ShB=Sh;
int[int] 	wall(3);     		wall    	= [1,7];			// 	fixed displacement

fespace VhSB2(ShB,P2),			VhSB1(ShB,P1),			VhSB0(ShB,P0)			;
VhSB0 		chi;		//	chi: characteristic function
VhSB0 		chiP;		//	chiP: characteristic function with weak domain

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

Mat A;
macro def(i)[i, i#B, i#C]// EOM     // vector field definition
macro init(i)[i, i, i]// EOM        // vector field initialization
createMat(Sh, A, Pk)
Wh def(u);
WhPlt [gu1, gu2, gu3];
WhPlt [pltx, plty, pltz];
WhPlt [reducex, reducey, reducez];
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,mu,0,0],[0,0,0,0,mu,0],[0,0,0,0,0,mu]] //elastic tensor
real gZ;	gZ = -0.01;
macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM
varf vPb(def(u), def(v)) = int3d(Sh)((D*epsilon(u))'*epsilon(v)*E*chiP)
						+int2d(Sh,3)(gZ*vC)
						+ on(wall, u = 0.0, uB = 0.0, uC = 0.0); // '

//-----------------------adding part------------------------
Mat B;
macro grad(phi)[dx(phi), dy(phi), dz(phi)]// EOM    // three-dimensional gradient
createMat(Sh, B, P1)
fespace WhRDE(Sh, P1);
WhRDE	phi,testphi;			// level set function and its test function
WhRDE	ophi;					// ophi: previous level set function
//-----------------------adding part------------------------

phi 		= 1.0;			//	initialize level set function
ophi 		= phi*(phi<=1.0)*(phi>=-1.0)+1.*(phi>1.0)-1.*(phi<-1.0); //  ophi is mapped to satisfy the upper and lower constraints.
plot(ophi,ShB,cmm="ophi",fill=1,wait=false,boundary=false,WindowIndex=0);
chi 	= (0.5+phi/matw*(15./16-phi^2/matw^2*(5./8-3./16*phi^2/matw^2)))*(phi>=-matw)*(phi<=matw)+1.*(phi>matw);
chiP  	= (1.-matd)*chi+matd;

real obj;				//objective function
real objPre;			//Previous objective function

for (int iter=0;iter<2;iter++) {
	
	A = vPb(Wh, Wh);
	real[int] rhs = vPb(0, Wh);
	set(A, sparams = "-pc_type gamg -ksp_view -ksp_max_it 200", bs = 3);
	u[] = A^-1 * rhs;
	
	if(!NoGraphicWindow) {
		real[int] sol;
		ChangeNumbering(A, u[], sol);
		ChangeNumbering(A, u[], sol, inverse = true);
		int[int] rest=restrict(Wh, WhPlt, n2o);
		for[i, v : rest] pltx[][v] = u[][i];
		mpiReduce(pltx[], reducex[], processor(0, mpiCommWorld), mpiSUM);
	}
	
	if(mpirank==0){
		gu1[]=reducex[];
		gu2[]=reducey[];
		gu3[]=reducez[];
		ShPltMoved = movemesh3(ShPlt,transfo=[x+gu1,y+gu2,z+gu3]);
		plot(ShPltMoved,cmm="Deformation");

		obj 		= int2d(ShPlt,3)(gZ*gu3);

		cout	<<"+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
		cout	<< "!!!	iter	ObjF "<<endl;
		cout	<< "!!!	"<<iter<<"	"<<obj<<endl;
		cout	<<"+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
		
		phi 	= 0.5*(cos((2*pi/1.0)*y)+0.3);			//	initialize level set function		
		cout	<<"+++++++++++++++++++++++++++++++++++++++++++++++++++"<=-1.0)+1.*(phi>1.0)-1.*(phi<-1.0); 
		chi 	= (0.5+phi/matw*(15./16-phi^2/matw^2*(5./8-3./16*phi^2/matw^2)))*(phi>=-matw)*(phi<=matw)+1.*(phi>matw);
		chiP  	=  (1.-matd)*chi+matd;
		plot(ophi,ShB,cmm="ophiO",fill=1,wait=false,boundary=false,WindowIndex=0);
		
}// End optimization loop

At this time, we obtained the following error.

Error line number 403, in file  macro: partitionPrivate in C:\Program Files (x86)\FreeFem++\\idp\macro_ddm.idp, before  token ;
Invalid array size  for  vectorial fespace function
  current line = 403 mpirank 0 / 8
Compile error : Invalid array size  for  vectorial fespace function
        line number :403, ;
error Compile error : Invalid array size  for  vectorial fespace function
        line number :403, ;
 code = 1 mpirank: 0

Is there any procedure required when two different PDEs coexist in the same code like this?

Sorry if there have been similar questions in the past.
I would appreciate your response.
Many thanks,
Keita Kambayashi

2

You need to scope the macro def and init definitions, cf. FreeFem-sources/stokes-block-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub.

3

Scoping a macro def? What language was this modeled from ?

4

C and C++ (among others).

5

Dear Professor Jolivet

Thank you very much for all your help.
First, I tried the scope. The calculation can be done by follow code:

verbosity = 3;
load "medit" 	// for data save
load "msh3"
load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"  
// Parameters
real 	E;			E 			= 1.0;						// Young's modulus	
real 	nu; 		nu 			= 0.3;						// Poisson's ratio
real 	lambda;		lambda 		= nu*1./(1+nu)/(1-2.*nu);	// Lame coefficient
real 	mu; 		mu 			= 1./2./(1+nu);				// Lame coefficient
real 	matd;		matd 		= 1e-3;
real 	matw;		matw		= 0.5;
real 	L;			L 			= 0.1;						// characteristic length
real 	tau; 		tau 		= 1.0e-6;					// regularization paramter
real 	CdF;		CdF 		= 0.7;						// normalization paramter for sensitivity
real 	dt;			dt			= 10;						// initial time increment
// Generating mesh
int 	paramesh;	paramesh	= 2;
int[int] 	labs(6);		// labels on boundaries of squre: (bottom, right, top, left)
mesh3 		Shc1;				// Mesh of the squre 1
labs		= [1,0,3,0,0,6];
Shc1 		= cube(5*paramesh,10*paramesh,1*paramesh,[-0.25+(0.25+0.25)*x, 0.0+(1.0-0.0)*y, -0.25+(-0.15+0.25)*z],label=labs,flags=1,region=1);
mesh3 		Shc2;				// Mesh of the squre  2
labs		= [7,8,9,100,11,12];	
Shc2 		= cube(5*paramesh,10*paramesh,4*paramesh,[-0.25+(0.25+0.25)*x, 0.0+(1.0-0.0)*y, -0.15+(0.25+0.15)*z],label=labs,flags=1,region=1);
mesh3 		Sh=Shc1+Shc2;//+Shc3;
mesh3		ShPlt=Sh;
mesh3 		ShPltMoved;
mesh3 		ShB=Sh;
int[int] 	wall(3);     		wall    	= [1,7];			// 	fixed displacement
int[int] 	phiOne(1);     		phiOne    	= [3];				//  phi=1.0
//define macros
macro 	defSTR(i)[i, i#B, i#C]	// EOM     	// vector field definition
macro 	initSTR(i)[i, i, i]		// EOM  	// vector field initialization
macro 	defRDE(i)i				// EOM     	// vector field definition
macro 	initRDE(i)i				// EOM     	// vector field initialization
func 	Pk = [P2, P2, P2]; // finite element space
func 	PkRDE = P1; // finite element space
// for restrict
int[int] n2o;
macro ShN2O()n2o// EOM
// fespaces
fespace Wh(Sh, Pk), WhPlt(ShPlt, Pk);
fespace WhRDE(Sh, PkRDE), WhPltRDE(ShPlt,P1);
// 
Mat A, B;
{
	buildDmesh(Sh);
	{
		macro def(u) defSTR(u)			// EOM     // vector field definition
		macro init(u) initSTR(u)		// EOM        // vector field initialization
		createMat(Sh, A, Pk)
	}
	{
		macro def(phi) defRDE(phi)	// EOM     // vector field definition
		macro init(phi) initRDE(phi)	// EOM        // vector field initialization
		createMat(Sh, B, PkRDE)
	}
}
Wh defSTR(u);
WhPlt [pltx, plty, pltz], [reducex, reducey, reducez], [gu1, gu2, gu3];
WhRDE defRDE(phi);
WhRDE		ophi, chi, chiP;		
WhPltRDE 	pltRDE, reduceRDE, gphi;
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,mu,0,0],[0,0,0,0,mu,0],[0,0,0,0,0,mu]] //elastic tensor
real 	gZ;		gZ = -0.01;
macro grad(phi)[ dx(phi), dy(phi), dz(phi) ]//EOM	// Three-dimensional gradient
// 3D elasticity problem
varf vPb(defSTR(u), defSTR(v)) = int3d(Sh)((D*epsilon(u))'*epsilon(v)*E*chiP)
						+int2d(Sh,3)(gZ*vC)+ on(wall, u = 0.0, uB = 0.0, uC = 0.0); // '
// Reaction Diffusion Equation
varf rde(defRDE(phi),defRDE(testphi))
				=int3d(Sh)((phi*testphi/dt)+( grad(phi)'*grad(testphi) )*tau*L^2)
				-int3d(Sh)((CdF*(10.0)+ophi/dt)*testphi)+on(phiOne,phi=1);
phi 		= 1.0;			//	initialize level set function
ophi 		= phi*(phi<=1.0)*(phi>=-1.0)+1.*(phi>1.0)-1.*(phi<-1.0); //  ophi is mapped to satisfy the upper and lower constraints.
chi 	= (0.5+phi/matw*(15./16-phi^2/matw^2*(5./8-3./16*phi^2/matw^2)))*(phi>=-matw)*(phi<=matw)+1.*(phi>matw);
chiP  	= (1.-matd)*chi+matd;
real obj;				//objective function
for (int iter=0;iter<3;iter++) {
		
		// solve 3D elasticity problem
		A = vPb(Wh, Wh);
		real[int] rhs = vPb(0, Wh);
		set(A, sparams = "-pc_type gamg -ksp_view -ksp_max_it 200", bs = 3);
		u[] = A^-1 * rhs;
		if(!NoGraphicWindow) {
			real[int] sol;
			ChangeNumbering(A, u[], sol);
			ChangeNumbering(A, u[], sol, inverse = true);
			int[int] rest=restrict(Wh, WhPlt, n2o);
			for[i, v : rest] pltx[][v] = u[][i];
			mpiReduce(pltx[], reducex[], processor(0, mpiCommWorld), mpiSUM);
		}
		if(mpirank==0){
			gu1[]=reducex[];
			gu2[]=reducey[];
			gu3[]=reducez[];
			ShPltMoved = movemesh3(ShPlt,transfo=[x+gu1,y+gu2,z+gu3]);
			plot(ShPltMoved,cmm="Deformation");
		}
		// solve Reaction Diffusion Equation
		B = rde(WhRDE, WhRDE);
		real[int] rhsRDE = rde(0, WhRDE);
		set(B, sparams = "-pc_type hypre -ksp_view -ksp_max_it 200");
		phi[] = B^-1 * rhsRDE;
		if(!NoGraphicWindow) {
			real[int] solRDE;
			ChangeNumbering(B, phi[], solRDE);
			ChangeNumbering(B, phi[], solRDE, inverse = true);
			int[int] restRDE=restrict(WhRDE, WhPltRDE, n2o);
			for[i, testphi : restRDE] pltRDE[][testphi] = phi[][i];
			mpiReduce(pltRDE[], reduceRDE[], processor(0, mpiCommWorld), mpiSUM);
		}
		if(mpirank==0){
			gphi[] = reduceRDE[];
			plot(gphi,ShPlt,cmm="global phi");
		}
		ophi 	= phi*(phi<=1.0)*(phi>=-1.0)+1.*(phi>1.0)-1.*(phi<-1.0); 
		chi 	= (0.5+phi/matw*(15./16-phi^2/matw^2*(5./8-3./16*phi^2/matw^2)))*(phi>=-matw)*(phi<=matw)+1.*(phi>matw);
		chiP  	=  (1.-matd)*chi+matd;			
		plot(ophi,ShPlt,cmm="ophi",fill=1,wait=false,boundary=false,WindowIndex=0);
}// End optimization loop

And here we have another important problem.

What we want to do:
On each processor, I’d like to solve a 3D elasticity problem based on the solution (chiP) of the RDE.

What is the issue:
The connection information between the divided regions by ddm (domain decomposition methods) is not good
(deformation is odd near the split face.).

How can this be avoided?
Do I need to save it once in a global solution and then distribute it?

Many thanks,
Kambayashi

6

The connection information between the divided regions by ddm (domain decomposition methods) is not good

What does that mean? Do you get a different solution with one process and multiple processes?

7

When solving a 3D elastic problem including RDE solution information,
in the overlapping surface of the divided area, an extremely small value appears.

8

Please provide a minimal working example showing such a small value. If this is only related to elasticity, there should be no need for the full code.

9

Dear Professor Jolivet

Thank you for your reply and patients.

I have modified the code to show minimal behavior.
But at that time, I have found cause that when the RDE is solved some times, the overlap surface produces values like numerical garbage.

verbosity = 3;
load "medit" 	// for data save
load "msh3"
load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"  
// Parameters
real 	matd;		matd 		= 1e-3;
real 	matw;		matw		= 0.5;
real 	L;			L 			= 0.1;						// characteristic length
real 	tau; 		tau 		= 1.0e-6;					// regularization paramter
real 	CdF;		CdF 		= 0.7;						// normalization paramter for sensitivity
real 	dt;			dt			= 10;						// initial time increment
// Generating mesh
int 	paramesh;	paramesh	= 2;
int[int] 	labs(6);		// labels on boundaries of squre: (bottom, right, top, left)
mesh3 		Shc1;				// Mesh of the squre 1
labs		= [1,0,3,0,0,6];
Shc1 		= cube(5*paramesh,10*paramesh,1*paramesh,[-0.25+(0.25+0.25)*x, 0.0+(1.0-0.0)*y, -0.25+(-0.15+0.25)*z],label=labs,flags=1,region=1);
mesh3 		Shc2;				// Mesh of the squre  2
labs		= [7,8,9,100,11,12];	
Shc2 		= cube(5*paramesh,10*paramesh,4*paramesh,[-0.25+(0.25+0.25)*x, 0.0+(1.0-0.0)*y, -0.15+(0.25+0.15)*z],label=labs,flags=1,region=1);
mesh3 		Sh=Shc1+Shc2;//+Shc3;
mesh3		ShPlt=Sh;
mesh3 		ShPltMoved;
mesh3 		ShB=Sh;
int[int] 	wall(3);     		wall    	= [1,7];			// 	fixed displacement
int[int] 	phiOne(1);     		phiOne    	= [3];				//  phi=1.0
//define macros
macro 	defRDE(i)i				// EOM     	// vector field definition
macro 	initRDE(i)i				// EOM     	// vector field initialization
func 	PkRDE = P1; // finite element space
// for restrict
int[int] n2o;
macro ShN2O()n2o// EOM
// fespaces
fespace WhRDE(Sh, PkRDE), WhPltRDE(ShPlt,P1);
// 
Mat B;
buildDmesh(Sh);
macro def(phi) defRDE(phi)	// EOM     // vector field definition
macro init(phi) initRDE(phi)	// EOM        // vector field initialization
createMat(Sh, B, PkRDE)

WhRDE defRDE(phi);
WhRDE		ophi, chi, chiP;		
WhPltRDE 	pltRDE, reduceRDE, gphi;
macro grad(phi)[ dx(phi), dy(phi), dz(phi) ]//EOM	// Three-dimensional gradient
// Reaction Diffusion Equation
varf rde(defRDE(phi),defRDE(testphi))
				=int3d(Sh)((phi*testphi/dt)+( grad(phi)'*grad(testphi) )*tau*L^2)
				-int3d(Sh)((CdF*(0.1*cos(2*pi*(y/0.5)))+ophi/dt)*testphi)+on(phiOne,phi=1);
phi 		= 1.0;			//	initialize level set function
ophi 		= phi*(phi<=1.0)*(phi>=-1.0)+1.*(phi>1.0)-1.*(phi<-1.0); //  ophi is mapped to satisfy the upper and lower constraints.
chi 	= (0.5+phi/matw*(15./16-phi^2/matw^2*(5./8-3./16*phi^2/matw^2)))*(phi>=-matw)*(phi<=matw)+1.*(phi>matw);
chiP  	= (1.-matd)*chi+matd;
for (int iter=0;iter<3;iter++) {
		// solve Reaction Diffusion Equation
		B = rde(WhRDE, WhRDE);
		real[int] rhsRDE = rde(0, WhRDE);
		set(B, sparams = "-pc_type hypre -ksp_view -ksp_max_it 200");
		phi[] = B^-1 * rhsRDE;
		if(!NoGraphicWindow) {
			real[int] solRDE;
			ChangeNumbering(B, phi[], solRDE);
			ChangeNumbering(B, phi[], solRDE, inverse = true);
			int[int] restRDE=restrict(WhRDE, WhPltRDE, n2o);
			for[i, testphi : restRDE] pltRDE[][testphi] = phi[][i];
			mpiReduce(pltRDE[], reduceRDE[], processor(0, mpiCommWorld), mpiSUM);
		}
		if(mpirank==0){
			gphi[] = reduceRDE[];
			plot(gphi,ShPlt,cmm="global phi");
		}
		ophi 	= phi*(phi<=1.0)*(phi>=-1.0)+1.*(phi>1.0)-1.*(phi<-1.0); 
		chi 	= (0.5+phi/matw*(15./16-phi^2/matw^2*(5./8-3./16*phi^2/matw^2)))*(phi>=-matw)*(phi<=matw)+1.*(phi>matw);
		chiP  	=  (1.-matd)*chi+matd;			
}// End optimization loop

This seems to be the cause of such small values.
How can this be resolved?

Many thanks,
Kambayashi

10

So if you solve this with an LU method you get the proper values?

11

Dear Professor Jolivet

Thank you for your reply.

I’m so sorry for mistake meaning.

iteratively → some times

12

Does the same problem appears with only a single process?

13

Dear Professor Jolivet

Thank you for your reply.
I’m so sorry for my delay response.

No, that problem doesn’t happen with only a single process.

Many thanks,
Kambayashi

14

You probably need to exchange() the ophi, chi, and chiP values to keep them in sync between processes.

15

Dear Professor Jolivet,

Thank you very much for all your replies.

While waiting for a response, I tried another method using 

 broadcast()

And now, I am investigating whether it is being executed correctly.
At the same time, I’d like to lean about using 
 exchange()

Now, it takes long time for checking my code, so my reply was delay.
I’m so sorry. I will report back to you when they are done.

Many thanks,
Kambayashi

16

Dear Professor. Jolivet

It’s been a while.

Now, I checking my code.
Then, I’d like to the meaning of the command exchange( ) you suggested.
I haven’t used it yet.

How can I know that meaning?
I usually search the code meaning in source file, but sometimes I can’t find or understand.
I’ll be very glad to your response about meaning and how to use exchange(), and how to find the meaning of the such codes for the future.

Thanks many in your busy schedules,
Keita Kambayashi

17
18

Dear Prof. Jolivet

Thank you for your reply!
I should check it in detail.

Best regards,
Keita Kambayashi