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
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);