Hello, I am experiencing issues with obtaining consistent results between serial and parallel computations of the three-dimensional heat conduction equation using PETSc. The first image is my mesh diagram, and the second is the temperature distribution diagram obtained after performing the calculation on a single core. The boundary conditions are the same. It can be observed that there is no heat generation effect inside the two cubes, which is incorrect. Could you please help me understand why this might be happening? I would greatly appreciate any advice you can provide. Thank you very much.
load “PETSc”
load “msh3”
load “medit”
load “tetgen”
load “iovtk”
macro dimension()3// EOM
macro Grad(u) [dx(u),dy(u),dz(u)] //EOM
macro INTERPLATE(x, x1, x2, y1, y2) ((x - x1) * (y1 - y2) / (x2 - x1) + y1) //
include “macro_ddm.idp”
real T0 = 300;
real voltage = 80;
real dt = 120;
string dir = “test”;
real localJoule = 0;
mesh3 Th = readmesh3(“test0.mesh”);
int[int]labs=labels(Th);
int[int]reg=regions(Th);
cout<<“lables\n”<<labs<<endl;
cout<<“regions\n”<<reg<<endl;
plot(Th, wait=true);
// define finite element
fespace W0(Th, P03d);
fespace W1(Th, P13d);
W1 Jex,Jey,Jez,Je2,Je,Joule;
W1 phi,phiv;
W1 T=T0,Tv=T0,Told=T0;
W0 kappa, rhoe, Cv;
real TE = int3d(Th, 3)(T) / int3d(Th, 3)(1);
real TH = int3d(Th, 2)(T) / int3d(Th, 2)(1);
real T1 = 300;
real T2 = 3800 + 300;
real b0kappa = 0.3;
real c0kappa = 400;
kappa = (region == 1)*b0kappa +
(region == 2)*c0kappa +
(region == 3)*c0kappa ;
real c0t3 = INTERPLATE(TH, T1, T2, 590, 100) / 2 * 1e-6;
real c8t25 = INTERPLATE(TE, T1, T2, 550, 100) / 2 * 1e-6;
real b0rhoe = 1e6;
rhoe = (region == 1)*b0rhoe +
(region == 2)*c0t3 +
(region == 3)*c8t25 ;
real c0cv = 590 * 1000;
real b0cv = 2000 * 858;
Cv = (region == 1)*b0cv +
(region == 2)*c0cv +
(region == 3)*c0cv ;
varf laplace(phi,phiv, eps=1.0e-8)=
int3d(Th)( Grad(phi)'*Grad(phiv)/rhoe)
+ on(2, phi=voltage/2)
+ on(3, phi=-voltage/2)
+ on(1, phi=0)
;
varf heat(T, Tv, eps=1.0e-8)=
int3d(Th)( Cv/dt * T * Tv )
-int3d(Th)( Cv/dt * Told * Tv )
+int3d(Th)( kappa * (Grad(T)'*Grad(Tv)) )
-int3d(Th)( Joule * Tv )
+on(1, T=T0);
Mat A ;
buildMat(Th, getARGV(“-split”, 1), A, P13d, mpiCommWorld)
real time01 = clock();
matrix Aloc= laplace(W1, W1);
A = Aloc;
real[int] b = laplace(0, W1);
real time02 = clock();
set(A, sparams = “-pc_type hypre -ksp_type gmres”);
phi = A^-1 * b;
real time03 = clock();
cout <<" Time_lap " << time02 - time01 << " s "
<<" Time_cal " << time03 - time02 << " s " << endl;
Jex = - dx(phi) / rhoe;
Jey = - dy(phi) / rhoe;
Jez = - dz(phi) / rhoe;
Je2 = pow(Jex,2) + pow(Jey,2) + pow(Jez,2);
Je = pow(Je2,0.5);
Joule = rhoe * Je2;
localJoule = int3d(Th)(Joule);
cout<<"Processor: “<<mpirank<<”:Local Joule Heat = "<<localJoule<<endl;
real globalJoule = 0.0;
mpiAllReduce(localJoule, globalJoule, mpiCommWorld, mpiSUM);
cout<<"Processor: “<<mpirank<<”:Total Joule Heat = "<<globalJoule<<endl;
// for (int i = 0; i < 11; i++){
Told = T;
real time04 = clock();
matrix Bloc= heat(W1, W1);
A = Bloc;
b = heat(0, W1);
real time05 = clock();
set(A, sparams = “-pc_type hypre -ksp_type gmres -ksp_rtol 1e-8”);
T = A^-1 * b;
real time06 = clock();
cout <<" Time_heat " << time05 - time04 << " s "
<<" Time_cal " << time06 - time05 << " s " << endl;
// }
string fn1 = dir+“/phi.txt”;
ofstream file(fn1);
file << phi << endl;
savevtk(dir+“/phi.vtk”, Th, phi, dataname=“phi”);
string fn2 = dir+“/T.txt”;
ofstream file1(fn2);
file1 << T << endl;
savevtk(dir+“/T.vtk”, Th, T, dataname=“T”);