The parallel solution of the 3D heat conduction equation using PETSc results in errors

yongheng · October 21, 2024, 2:20pm

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

prj · October 21, 2024, 2:59pm

I can’t copy/paste your code, so I can’t help you.

prj · October 21, 2024, 3:02pm

I see at least two obvious issues:

you are using the default tgv value of 1e+30, which is bad for iterative solvers;
you are integrating a function on Th which has overlapping elements, and then summing all contributions without taking into account multiplicities, so the value can’t be right.

yongheng · October 21, 2024, 3:23pm

Thank you very much for your answer. I tried setting tgv=-1, but got the same incorrect result. There is a strange phenomenon that the values obtained from most nodes are opposite in sign to the serial values.

prj · October 21, 2024, 3:28pm

Do you get the same results with the PETSc code and only a single MPI process?

fb77 · October 21, 2024, 4:24pm

I think that you are falling into the trap of varf-on, see

For your heat definition you should probably change the sign of the rhs parts:
+int3d(Th)( Cv/dt * Told * Tv )
+int3d(Th)( Joule * Tv )

yongheng · October 22, 2024, 3:09am

I apologize, but when I used fb77’s approach, the results obtained with a single MPI process were consistent, but I still couldn’t get the correct results when using multiple MPI processes. Here is my mesh file.
test.mesh (1.9 MB)
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 +kappa * (Grad(T)'*Grad(Tv)) )
+int3d(Th)( Cv/dt * Told * Tv + 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, tgv = -1);
A = Bloc;
b = heat(0, W1, tgv = -1);
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”);

prj · October 22, 2024, 4:08am

You are still not addressing any of the two issues I raised earlier. I still can’t copy/paste your code.

yongheng · October 22, 2024, 9:16am

Here are my code file and mesh file. I have set tgv to -1, but I still don’t know how to modify the second issue you mentioned. Could you take a look at it for me?
test0.edp (3.8 KB)
test.mesh (1.9 MB)

prj · October 22, 2024, 9:39am

You can have a look here FreeFem-sources/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub on how to integrate a function on a non-overlapping mesh to get the proper results in parallel.

Topic		Replies	Views
Solve 3D thermoelasticity with PETSc General Discussion	6	163	May 24, 2024
Poor results with PETSc General Discussion	4	430	June 24, 2021
Computing numerical solution with parallel version General Discussion	3	303	March 4, 2022
Unsteady three dimensional PETSc General Discussion	6	773	May 27, 2021
PETSc on 3D Cube with Periodic Boundary Conditions General Discussion	6	220	November 14, 2023

The parallel solution of the 3D heat conduction equation using PETSc results in errors

Related topics