Solve 3D thermoelasticity with PETSc

Hello there !

I am trying to solve a 3D thermoelastic system on a cylinder but I encounter a few problems with the use of PETSc (and perhaps iovtk).

Indeed, I do not obtain the same mechanical results when I use more processors than one (thermic seems fine however).

In my program, I need to do some postprocessing on the mechanical displacement that requires some summation on the global mesh and I also need to have global variables for the different functional variables (Temperature and displacement) and stress tensor components.

What are the problems in my program ?

I put a copy of the code below.

Thank you in advance for your help.


load "msh3"
load "iovtk"
load "medit"
load "metis"
load "MUMPS"
load "mmg"
load "hpddm"
load "PETSc"
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions
macro def2(i)[i, i#B, i#C]//
macro init2(i)[i, i, i]//

verbosity =1;
int master = 0;

mesh3 Th,Sh,Thlocal,Shlocal ;
int Robin=3; int Neumann=2; int fixe=1; int interface=0;int fixebas=4;

int N=50;//50
border cc(t=0,2*pi){x=1.*cos(t);y=1.*sin(t);label=Robin;}
mesh Th2= buildmesh(cc(N));
real z0=0. ;
real z1 = 1.;
int[int] labelmidvec = [0,Robin,1,Robin,2,Robin,3,Robin,4,Robin];
int[int] labelupvec = [0,fixe,1,fixe,2,fixe,3,fixe,4,fixe];
int[int] labeldownvec = [0,fixebas,1,fixebas,2,fixebas,3,fixebas,4,fixebas];

Th = buildlayers(Th2, N, zbound=[z0,z1],
  labelmid=labelmidvec, labelup = labelupvec, labeldown = labeldownvec);
func Pk = P1;
func Pkm = [P1, P1, P1];
int[int] n2o;
macro ThN2O()n2o// EOM

/* Local FE Spaces */
fespace ScalarField0(Th,P0);
fespace ScalarField1(Th,Pk);
fespace ScalarField2(Th, Pk);
fespace VectorField1(Th,Pkm);
fespace VectorField2(Th,Pkm);

/* Global FE Spaces */
fespace ScalarFieldSh0(Sh,P0);
fespace ScalarFieldSh1(Sh,Pk);
fespace VectorFieldSh1(Sh,Pkm);

/*From local to Global*/
int[int] R0 = restrict(ScalarField0, ScalarFieldSh0, n2o);
int[int] R1 = restrict(ScalarField1, ScalarFieldSh1, n2o);
int[int] RV1 = restrict(VectorField1, VectorFieldSh1, n2o);
ScalarFieldSh0 hrtemp ;
ScalarFieldSh1 Ttemp ;
VectorFieldSh1 def2(utemp) ;

/*Matrix PETSc*/

Mat A0;
createMat(Th, A0, P0);

Mat A;
createMat(Th, A, Pk);

/*Mechanical matrix*/

Mat Am;
macro def(u)def2(u)//
macro init(u)init2(u)//
createMat(Th, Am, Pkm);

ScalarField1 T; // Temperature
ScalarField1 Trenorm; // renormalized temperature, wrt the boundary
ScalarField1 Tst; // adjoint temperature
ScalarField1 vSh; //adjoint temperature test

VectorField1 def2(u); // Displacement
VectorField1 def2(ust); // adjoint displacement
VectorField1 def2(utilde); // displacement test function
ScalarField1 eps11, eps22, eps33, eps13,eps12,eps23; // deformation
ScalarField0 sigma11,sigma13,sigma12,sigma23,sigma22,sigma33; // stress
ScalarField0 sigmarr, sigmatt;	// radial and angular stress

ScalarFieldSh0 sigmaGlob11,sigmaGlob13,sigmaGlob12,sigmaGlob23,sigmaGlob22,sigmaGlob33;
ScalarFieldSh1 TrenormGlob,TGlob,TstGlob;
VectorFieldSh1 def2(uGlob),def2(ustGlob);

// Operators
real sqrt2 = sqrt(2.);

// gradient 3D of scalar function
macro grad3D(u)[dx(u), dy(u),dz(u)] // EOM

// symmetric gradient 3D of vector field
macro epsilon3D(u) [dx(u), dy(u#B), dz(u#C),(dz(u) + dx(u#C)),(dz(u#B) + dy(u#C)), (dy(u) + dx(u#B))]  // EOM
// symmetric gradient 3D of vector field
macro epsilon3Dloc(u) [dx(u), dy(u#B), dz(u#C),(dz(u) + dx(u#C)),(dz(u#B) + dy(u#C)), (dy(u) + dx(u#B))] // EOM

// divergence 3D of vector field
macro div3D(u) (dx(u) + dy(u#B) + dz(u#C)) // EOM

// divergence 3D of vector field
macro div3Dloc(u) (dx(u) + dy(u#B) + dz(u#C)) // EOM

func real sigma11F(real eps11, real eps22, real eps33, real eps13, real eps23,real eps12, real T){
	return (2.*mu+lambda)*eps11+lambda*eps22+lambda*eps33-(2.*mu+3.*lambda)*alpha*(T - TambiantPhys);

///	PDE	///

/*Thermal PDE*/
varf heatvarf(T,v) = 
		int3d(Th)((1. * grad3D(T))' * grad3D(v))
	    + int3d(Th)(1. * v) 
		+ int2d(Th,3)(1.*T*v)
		+ int2d(Th,3)(1.*Tborder*v); 

A = heatvarf(ScalarField1, ScalarField1);
set(A, sparams = "-pc_type hypre -ksp_type gmres -ksp_max_it 200");
real[int] bheat = heatvarf(0, ScalarField1);
Trenorm[] = A^-1 * bheat;

	real[int] tmpT;	
	ChangeNumbering(A, Trenorm[], tmpT);
	ChangeNumbering(A, Trenorm[], tmpT, inverse = true);
	Ttemp[](R1) = Trenorm[];
	mpiReduce(Ttemp[], TrenormGlob[], processor(0, mpiCommWorld), mpiSUM);
	TGlob = TrenormGlob;

/*Mechanical problem*/

varf elasticitevarf(def2(u),def2(utilde)) =
   int3d(Th)((2.*1. * epsilon3D(u))' * epsilon3D(utilde))
	+ int3d(Th)((1. * div3D(u))' * div3D(utilde))
		 + int3d(Th)(1.e-10*1.*(u'*utilde))
	+int3d(Th)((T-1.)*(2.*1.+3.*1.)*1.*(div3D(utilde)) )
	+ on(fixebas,uC=0);

Am = elasticitevarf(VectorField1, VectorField1);
set(Am, sparams = "-pc_type lu -pc_factor_mat_solver_type mumps");
real[int] bmeca = elasticitevarf(0, VectorField1);
u[] = Am^-1 * bmeca;

def2(uGlob)=def2(uGlob) - def2(uGlob);

	real[int] tmpu;
	ChangeNumbering(Am, u[], tmpu);
	ChangeNumbering(Am,u[], tmpu, inverse = true);
	utemp[](RV1) = u[];
	mpiReduce(utemp[], uGlob[], processor(0, mpiCommWorld), mpiSUM);

real moyenneu1 = int3d(Sh)(uGlob) / 1.;
real moyenneu2 = int3d(Sh)(uGlobB) / 1.;
real moyennediffderivee= int3d(Sh)(dy(uGlob)-dx(uGlobB)) / 1.;
VectorFieldSh1 def2(utest) = [uGlob - moyenneu1  - 0.5 * moyennediffderivee * y, 
		uGlobB - moyenneu2 +  0.5 * moyennediffderivee * x, 
u[] = uGlob[](RV1);

eps11 = dx(u);
eps22 = dy(uB);
eps33 = dz(uC);
eps12 =  (dy(u) + dx(uB))/2.;
eps13 =  (dz(u) + dx(uC))/2.;
eps23 =  (dz(uB) + dy(uC))/2.;

sigma11 = sigma11F(eps11, eps22, eps33, eps13, eps23, eps12, T);

sigmaGlob11 =0. ;

real[int] tmps11;	
ChangeNumbering(A0, sigma11[], tmps11);
ChangeNumbering(A0, sigma11[], tmps11, inverse = true);
hrtemp[](R0) = sigma11[];	
mpiReduce(hrtemp[], sigmaGlob11[], processor(0, mpiCommWorld), mpiSUM);

ScalarFieldSh1 Tdisp = TGlob;
int[int] order=[1]; 
savevtk("Temperature.vtu", Sh,  Tdisp, dataname=" Temperature", order=order, append=true);

ScalarFieldSh1 uDisp =uGlob;
ScalarFieldSh1 uBDisp =uGlobB;
ScalarFieldSh1 uCDisp =uGlobC;
savevtk("u.vtu", Sh,  uDisp, dataname=" u", order=order, append=true);
savevtk("uB.vtu", Sh,  uBDisp, dataname=" uB", order=order, append=true);
savevtk("uC.vtu", Sh,  uCDisp, dataname=" uC", order=order, append=true);

ScalarFieldSh1 sigma11Disp =sigmaGlob11;
savevtk("sigma11.vtu", Sh, sigma11Disp, dataname="sigma11", order=order, append=true);

What are the problems in my program ?

That is for you to tell us.

My problem is that when I use more than one processor, I obtain false results for the displacement field u and the stress sigma11 but I do not know why it is the case.

I think that it could come from a poor resolution of the elasticity equation or from the transition of the local displacement to the global displacement. But I am not certain and I do not see the problem in my code.

OK, the first issue is that I simply cannot run your code. Please send a runnable snippet.

Here is the runnable file :
LocalGlobal.edp (6.8 KB)
I will also try to modify the original post to be runnable.

That file is 250-line long. If your problem comes from a linear solve, it should only be around 20 lines of code. Please send a minimal working example reproducing the problem.

I will prepare a smaller code and try to find more precisely the location of the problem(s) in my code. Thank you for your time