# Solve 3D thermoelasticity with PETSc
**URL:** https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264
**Category:** General Discussion
**Created:** [May 23, 2024, 7:16am UTC](https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264 "2024-05-23T07:16:51Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Timothee99](https://community.freefem.org/u/Timothee99)
#### Post date: [May 23, 2024, 7:16am UTC](https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264/1 "2024-05-23T07:16:51Z")
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.
Timothée
```auto
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);
Sh=Th;
func Pk = P1;
func Pkm = [P1, P1, P1];
int[int] n2o;
macro ThN2O()n2o// EOM
buildDmesh(Th);
/* 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;
T=Trenorm;
TrenormGlob=0.;
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,
uGlobC];
def2(uGlob)=def2(utest);
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);
```
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [May 24, 2024, 5:17am UTC](https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264/2 "2024-05-24T05:17:52Z")
> What are the problems in my program ?
That is for you to tell us.
---
### Author:  [@Timothee99](https://community.freefem.org/u/Timothee99)
#### Post date: [May 24, 2024, 6:26am UTC](https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264/3 "2024-05-24T06:26:23Z")
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.
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [May 24, 2024, 6:39am UTC](https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264/4 "2024-05-24T06:39:05Z")
OK, the first issue is that I simply cannot run your code. Please send a runnable snippet.
---
### Author:  [@Timothee99](https://community.freefem.org/u/Timothee99)
#### Post date: [May 24, 2024, 6:56am UTC](https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264/5 "2024-05-24T06:56:32Z")
Here is the runnable file :
[LocalGlobal.edp](https://community.freefem.org/uploads/short-url/jZrVB08ig6D6ilvaOmKOPBwaeFq.edp) (6.8 KB)
I will also try to modify the original post to be runnable.
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [May 24, 2024, 9:31am UTC](https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264/6 "2024-05-24T09:31:41Z")
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.
---
### Author:  [@Timothee99](https://community.freefem.org/u/Timothee99)
#### Post date: [May 24, 2024, 3:57pm UTC](https://community.freefem.org/t/solve-3d-thermoelasticity-with-petsc/3264/7 "2024-05-24T15:57:21Z")
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