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

**URL:** <https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555>\
**Category:** General Discussion\
**Created:** [October 21, 2024, 2:20pm UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555 "2024-10-21T14:20:10Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![yongheng](https://avatars.discourse-cdn.com/v4/letter/y/f04885/32.png) [@yongheng](https://community.freefem.org/u/yongheng)\
**Post date:** [October 21, 2024, 2:20pm UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/1 "2024-10-21T14:20:10Z")

</div>

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.

 ![Screenshot from 2024-10-21 22-08-30](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/d/db7d6737eb138295ec02e529dd4e2dfbd5dfe0b8.png)  
 ![Screenshot from 2024-10-21 22-09-53](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/d/d6d28f03b95b4be7a2122548de695031439f2e3c.jpeg)  
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”);

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [October 21, 2024, 2:59pm UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/2 "2024-10-21T14:59:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [October 21, 2024, 3:02pm UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/3 "2024-10-21T15:02:25Z")

</div>

I see at least two obvious issues:

1. you are using the default `tgv` value of 1e+30, which is bad for iterative solvers;
2. 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.

---

<div class="post-metadata">

**Author:** ![yongheng](https://avatars.discourse-cdn.com/v4/letter/y/f04885/32.png) [@yongheng](https://community.freefem.org/u/yongheng)\
**Post date:** [October 21, 2024, 3:23pm UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/4 "2024-10-21T15:23:21Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [October 21, 2024, 3:28pm UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/5 "2024-10-21T15:28:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [October 21, 2024, 4:24pm UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/6 "2024-10-21T16:24:02Z")

</div>

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

> [@Change of sign in linear part of variational problem (varf vs problem)](https://community.freefem.org/t/change-of-sign-in-linear-part-of-variational-problem-varf-vs-problem/2984/2):
>
> solve vSolve(u, v) = a(u, v) + l(v) solves a(u, v) + l(v) = 0 while varf vPb(u, v) = a(u, v) + l(v) solves a(u, v) = l(v), so it is expected that you need to shift the sign in front of the linear form.

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

---

<div class="post-metadata">

**Author:** ![yongheng](https://avatars.discourse-cdn.com/v4/letter/y/f04885/32.png) [@yongheng](https://community.freefem.org/u/yongheng)\
**Post date:** [October 22, 2024, 3:09am UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/8 "2024-10-22T03:09:28Z")

</div>

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](https://community.freefem.org/uploads/short-url/54utbBABKdfXJ0CWUDPmzEQwH5g.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”);

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [October 22, 2024, 4:08am UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/9 "2024-10-22T04:08:33Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![yongheng](https://avatars.discourse-cdn.com/v4/letter/y/f04885/32.png) [@yongheng](https://community.freefem.org/u/yongheng)\
**Post date:** [October 22, 2024, 9:16am UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/10 "2024-10-22T09:16:51Z")

</div>

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](https://community.freefem.org/uploads/short-url/exqxaKC2McrpyE4HD705QxNaNCN.edp) (3.8 KB)  
[test.mesh](https://community.freefem.org/uploads/short-url/54utbBABKdfXJ0CWUDPmzEQwH5g.mesh) (1.9 MB)

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [October 22, 2024, 9:39am UTC](https://community.freefem.org/t/the-parallel-solution-of-the-3d-heat-conduction-equation-using-petsc-results-in-errors/3555/11 "2024-10-22T09:39:27Z")

</div>

You can have a look here [FreeFem-sources/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp#L17-L30) on how to integrate a function on a non-overlapping mesh to get the proper results in parallel.