# MOLFiniteDifference for 2D transient Heat equation with heat generation

**URL:** <https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498>\
**Category:** General Usage\
**Tags:** question\
**Created:** [January 7, 2025, 1:21pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498 "2025-01-07T13:21:59Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Reza\_AO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reza_ao/32/51137_2.png) [@Reza\_AO](https://discourse.julialang.org/u/Reza_AO)\
**Post date:** [January 7, 2025, 1:21pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/1 "2025-01-07T13:21:59Z")

</div>

Hi  
I am trying to solve the underlying **system** :

```julia
eq = [ρ*Cₚ*Dt(T(x,y,t)) ~ kₓ*Dxx(T(x,y,t)) + ky*Dyy(T(x,y,t)) + H*W*Re ,
      Re ~ A*exp(-E/(R*T(x,y,t)))*c(x,y,t)]

```

The **Domain** is :

```julia
domains = [x ∈ Interval(x_min, 0.013),
        y ∈ Interval(0,0.065),
        t ∈ Interval(0, 60*100)]

```

the system with the **parameters** :

 ![1](https://global.discourse-cdn.com/julialang/original/3X/f/9/f9462dc230be47fc91102c58c26cbea5f02e1018.jpeg)

The initial and Bondary conditions are:

 ![2](https://global.discourse-cdn.com/julialang/original/3X/8/9/89fa897808c716878b79818f0b820bea3ccbd47a.jpeg)

**Discretization** an **solve** :

```julia
N = 40
M = 20
order = 2
discretization = MOLFiniteDifference([x=>N, y=>M], t, approx_order=order)
@time prob = discretize(pdesys,discretization)
@time sol = solve(prob, AutoTsit5(Rosenbrock23()), saveat = 1)

```

When I solve this problem either **scaled or not** , The **shape of the mean** values of my outputs are correct but they have a delay as shown below:

 ![4](https://global.discourse-cdn.com/julialang/original/3X/e/4/e4bf4518ddb8885563907c5cb4070468c4ff80e7.jpeg)

the delay also exists for **c(x,y,t)**.  
I have used **small** N,M and **small time steps** but the results are the same.  
what am I doing wrong in here?  
PS: Because there are Neumann and Robin BCs in this problem, **The corner points of my problem resulted 0** ; Does this make sense? Could these 0 corners be the main reason for the delay?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 7, 2025, 2:36pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/2 "2025-01-07T14:36:13Z")

</div>

Can you share the full code for recreating that graph with both forms?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [January 7, 2025, 5:32pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/3 "2025-01-07T17:32:47Z")

</div>

You have such widely different parameter values. You can probably rescale t,x and y to reduce it, no?

---

<div class="post-metadata">

**Author:** ![Reza\_AO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reza_ao/32/51137_2.png) [@Reza\_AO](https://discourse.julialang.org/u/Reza_AO)\
**Post date:** [January 7, 2025, 9:53pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/4 "2025-01-07T21:53:30Z")

</div>

Thanks for your reply.  
tim and T\_exact are imported from an Excel file.

Code with dimensions:

```julia
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets,Statistics
using XLSX
using Plots
using Plots.PlotMeasures
custom_font = "Times bold"

@parameters x y t
@variables T(..) c(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

x_max = 0.013 #m
y_max = 0.065 #m
ρ = 2231.2 #Density
Cₚ = 1100 #[J/(kg*K)]Heat capacity
kₓ = 0.7 #[W/(m*K)] Thermal conductivity in x dir
ky = 140 #[W/(m*K)] Thermal conductivity in y dir
h_conv = 20 #[W/(m^2*K)] heat transfer oefficient
Tᵢ = 20 + 273.15 #[K] Initial temperature
Tₐ = 250 + 273.15 #[K]
αₓ = kₓ/(ρ*Cₚ) 
αy = ky/(ρ*Cₚ)
A = 5.14e25 #[1/s]	5.14E25 1/s
E = 2.7e5 #[J/mol]	2.7E5 J/mol
H = 6.2e5 #[J/kg]	6.2E5 J/kg
W = 334.68 #[kg/m^3]	334.68 kg/m³
R = 8.314
t_max = 100*60 
Bi = (h_conv*x_max)/kₓ
Re = -Dt(c(x,y,t))

eq = [ρ*Cₚ*Dt(T(x,y,t)) ~ kₓ*Dxx(T(x,y,t)) + ky*Dyy(T(x,y,t)) + H*W*Re ,
      Re ~ A*exp(-E/(R*T(x,y,t)))*c(x,y,t)]

domains = [x ∈ Interval(0, x_max),
        y ∈ Interval(0,y_max),
        t ∈ Interval(0, t_max)]

bcs = [T(x,y,0) ~ Tᵢ,
      Dx(T(0,y,t)) ~ 0, 
      kₓ*Dx(T(x_max,y,t)) ~ h_conv * (Tₐ - T(x_max,y,t)),
      Dy(T(x,0,t)) ~ 0, 
      Dy(T(x,y_max,t)) ~ 0,
      c(x,y,0) ~ 1,
      Dt(c(x,y,0)) ~ 0,
      Dx(c(0,y,t)) ~ 0,
      Dx(c(x_max,y,t)) ~ 0,
      Dy(c(x,0,t)) ~ 0,
      Dy(c(x,y_max,t)) ~ 0
]

@named pdesys = PDESystem(eq, bcs, domains,[x,y,t],[T(x, y,t),c(x,y,t)])

N = 40
M = 20
order = 2
discretization = MOLFiniteDifference([x=>N, y=>M], t, approx_order=order)
@time prob = discretize(pdesys,discretization)
@time sol = solve(prob, AutoTsit5(Rosenbrock23()), saveat = 1)

discrete_x = sol[x]
discrete_y = sol[y]
discrete_t = sol[t]

solT = sol[T(x, y, t)]
solc = sol[c(x, y, t)]

aveT = mean(solT,dims=[1,2])
aveC = mean(solc,dims=[1,2])

using Plots.PlotMeasures
custom_font = "Times bold"

plot(tim,T_exact , linewidth=3, alpha=1, color=:red,
 linestyle=:solid, xlabel="Time (min)",
  ylabel="Temperature (K)", guidefontsize=14,tickfontsize=14, legendfontsize=14,label="Real",
  guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
  left_margin = [6mm 0mm],
  bottom_margin = [6mm 0mm],
  dpi=300)
plot!(discrete_t/60, aveT[1,1,:] , linewidth=3, alpha=1, color=:black,
 linestyle=:solid, xlabel="Time (min)",
  ylabel="Temperature (K)", guidefontsize=14,tickfontsize=14, legendfontsize=14,label="MOL",
  guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
  left_margin = [6mm 0mm],
  bottom_margin = [6mm 0mm],
  dpi=300)

```

Result:

 ![11](https://global.discourse-cdn.com/julialang/original/3X/0/7/07d043165a66d0ede29855c480f2f79f692109bb.jpeg)

Code with scaling:

```julia
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets,Statistics
using XLSX
using Plots
using Plots.PlotMeasures
custom_font = "Times bold"

@parameters x y t
@variables T(..) c(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

x_max = 0.013 # width m
y_max = 0.065 # length m 
ρ = 2231.2 #Density
Cₚ = 1100 #[J/(kg*K)]Heat capacity
kₓ = 0.7 #[W/(m*K)] Thermal conductivity in x dir
ky = 140 #[W/(m*K)] Thermal conductivity in y dir
h_conv = 20 #[W/(m^2*K)] heat transfer oefficient
Tᵢ = 20 + 273.15 #[K] Initial temperature
Tₐ = 250 + 273.15 #[K]	Ambient temperature
αₓ = kₓ/(ρ*Cₚ) 
αy = ky/(ρ*Cₚ)
A = 5.14e25 #[1/s]	5.14E25 1/s
E = 2.7e5 #[J/mol]	2.7E5 J/mol
H = 6.2e5 #[J/kg]	6.2E5 J/kg
W = 334.68 #[kg/m^3]	334.68 kg/m³
R = 8.314
T_max = 370.72 + 273.15 # characteristic Temperature: maximum temperature from real solution
t_max = 27.835*60 # characteristic time: The time that maximum temperature occure
ΔT = T_max - Tᵢ
θₐ = (Tₐ - Tᵢ) /ΔT # scaled ambient temperature
θᵢ = 0 # scaled initial temperature
tf = 100*60/t_max # scaled final time
Re = -Dt(c(x,y,t))

eq = [Dt(T(x,y,t)) ~ (αₓ/(x_max^2))*t_max*Dxx(T(x,y,t)) + (αy/(y_max^2))*t_max*Dyy(T(x,y,t)) + (H*W*Re)/(ρ*Cₚ*ΔT) ,
      Re ~ A*t_max*exp(-E/(R*(T(x,y,t)*ΔT+Tᵢ)))*c(x,y,t)]

domains = [x ∈ Interval(0, 1),
        y ∈ Interval(0,1),
        t ∈ Interval(0, tf)]

bcs = [T(x,y,0) ~ θᵢ,
      Dx(T(0,y,t)) ~ 0, 
      Dx(T(1,y,t)) ~ ((h_conv*x_max)/kₓ) * (θₐ - T(1,y,t)),
      Dy(T(x,0,t)) ~ 0, 
      Dy(T(x,1,t)) ~ 0,
      c(x,y,0) ~ 1,
      Dt(c(x,y,0)) ~ 0,
      Dx(c(0,y,t)) ~ 0,
      Dx(c(1,y,t)) ~ 0,
      Dy(c(x,0,t)) ~ 0,
      Dy(c(x,1,t)) ~ 0
]

@named pdesys = PDESystem(eq, bcs, domains,[x,y,t],[T(x, y,t),c(x,y,t)])

N = 40
M = 20
order = 2
discretization = MOLFiniteDifference([x=>N, y=>M], t, approx_order=order)
@time prob = discretize(pdesys,discretization)
@time sol = solve(prob, AutoTsit5(Rosenbrock23()), saveat = 0.01)

discrete_x = sol[x]
discrete_y = sol[y]
discrete_t = sol[t]

solT = sol[T(x, y, t)]
solc = sol[c(x, y, t)]

aveT = mean(solT,dims=[1,2])
aveC = mean(solc,dims=[1,2])

plot(tim,T_exact , linewidth=3, alpha=1, color=:red,
 linestyle=:solid, xlabel="Time (min)",
  ylabel="Temperature [K]", guidefontsize=14,tickfontsize=14, legendfontsize=14,label="Real",
  guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
  left_margin = [6mm 0mm],
  bottom_margin = [6mm 0mm],
  dpi=300)

temp = aveT[1,1,:] * ΔT .+ Tᵢ
time_ = discrete_t*t_max/60
plot!(time_, temp , linewidth=3, alpha=1, color=:black,
 linestyle=:solid, xlabel="Time (min)",
  ylabel="Temperature (K)", guidefontsize=14,tickfontsize=14, legendfontsize=14,label="MOL",
  guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
  left_margin = [6mm 0mm],
  bottom_margin = [6mm 0mm],
  dpi=300)

```

Results:

 ![33](https://global.discourse-cdn.com/julialang/original/3X/0/4/0415cf6faff5645c3ada34ee6e7bddd129301c17.jpeg)

 ![22](https://global.discourse-cdn.com/julialang/original/3X/7/8/78f982509d2f304a7d86aa14bd21f1cc7f0c72f8.jpeg)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 8, 2025, 7:17am UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/6 "2025-01-08T07:17:47Z")

</div>

Can you share the code for how the values from the excel sheet are produced?

---

<div class="post-metadata">

**Author:** ![Reza\_AO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reza_ao/32/51137_2.png) [@Reza\_AO](https://discourse.julialang.org/u/Reza_AO)\
**Post date:** [January 8, 2025, 8:08am UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/7 "2025-01-08T08:08:37Z")

</div>

So I have 2 Excel files:  
one is the mean temperature and the other is the maximum temperature. I used the Maximum file just for the scaling using T\_max and t\_max. For mean temperature I just plotted the results as below.  
The mean Temperature file:

[https://docs.google.com/spreadsheets/d/1Tayb-kyOgOtMbRe1cvUWie37lPa3oMNL/edit?usp=sharing&ouid=100996059526622735141&rtpof=true&sd=true](https://docs.google.com/spreadsheets/d/1Tayb-kyOgOtMbRe1cvUWie37lPa3oMNL/edit?usp=sharing&ouid=100996059526622735141&rtpof=true&sd=true)

```julia
# Load and process data from the Excel file
xc = XLSX.readxlsx("avetemp.xlsx")
sh1 = xc["Sheet1"]
tim = sh1[1,4:end]
tim = vec(tim)
tim = convert(Vector{Float64}, tim)
T_exact = sh1[2,4:end]
T_exact = vec(T_exact)
T_exact = convert(Vector{Float64}, T_exact)
plot(tim,T_exact , linewidth=3, alpha=1, color=:red,label = "Mean Temperature",
 linestyle=:solid, xlabel="Time (min)",
  ylabel="Temperature (K)", guidefontsize=14,tickfontsize=14, legendfontsize=14,
  guidefont=custom_font, tickfont=custom_font, legendfont=custom_font,framestyle=:box,
  left_margin = [6mm 0mm],
  bottom_margin = [6mm 0mm],
  dpi=300)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 8, 2025, 8:54am UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/8 "2025-01-08T08:54:32Z")

</div>

Yes but can you share what you’re using as the ground truth for the datagen there?

---

<div class="post-metadata">

**Author:** ![Reza\_AO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reza_ao/32/51137_2.png) [@Reza\_AO](https://discourse.julialang.org/u/Reza_AO)\
**Post date:** [January 8, 2025, 5:13pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/9 "2025-01-08T17:13:19Z")

</div>

It’s from COMSOL multiphysics.

The real problem is from the paper below:  
[https://www.sciencedirect.com/science/article/abs/pii/S037877531830819X](https://www.sciencedirect.com/science/article/abs/pii/S037877531830819X)

I have validated full model via COMSOL. then For simplicity, I just considered the electrolyte to generate heat.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 9, 2025, 8:30am UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/10 "2025-01-09T08:30:41Z")

</div>

How do you know it’s not a property of the discretization?

---

<div class="post-metadata">

**Author:** ![Reza\_AO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reza_ao/32/51137_2.png) [@Reza\_AO](https://discourse.julialang.org/u/Reza_AO)\
**Post date:** [January 9, 2025, 1:04pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/11 "2025-01-09T13:04:06Z")

</div>

My apologies  
I’m not able to understand what you mean.  
Is there a problem with MOLFD discretization?  
COMSOL uses FEM (shape functions and orders).

---

<div class="post-metadata">

**Author:** ![Reza\_AO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reza_ao/32/51137_2.png) [@Reza\_AO](https://discourse.julialang.org/u/Reza_AO)\
**Post date:** [January 10, 2025, 2:09pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/12 "2025-01-10T14:09:46Z")

</div>

Could the cause of the problem be the high exponential term that changes rapidly?  
the heat generation caused by Re is highly nonlinear (from 0, it goes up to 4.5e5 in 4 minutes then falls to 0 again). Must I choose Nonelinear solvers?

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [January 10, 2025, 2:31pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/13 "2025-01-10T14:31:16Z")

</div>

> [@Reza\_AO](#):
>
> **The corner points of my problem resulted 0** ; Does this make sense? Could these 0 corners be the main reason for the delay?

I think these points are never used in the discretization and cannot be responsible for the delay. But they could be responsible for a small amplitude error when you calculate/diagnose the mean values later on. You should calculate the mean values only on the interior points.

BTW, I notice that nothing in your problem depends on y. Maybe the example provided here is a simplified version of a more complicated setup with an explicit y dependence, otherwise you could work in 1D I suppose?

---

<div class="post-metadata">

**Author:** ![Reza\_AO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reza_ao/32/51137_2.png) [@Reza\_AO](https://discourse.julialang.org/u/Reza_AO)\
**Post date:** [January 10, 2025, 4:59pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/14 "2025-01-10T16:59:42Z")

</div>

Thanks for Your reply. Yes it’s simplified. I don’t think the problem is from those zero points because even the maximum temperature has the delay.

---

<div class="post-metadata">

**Author:** ![Stephen](https://avatars.discourse-cdn.com/v4/letter/s/96bed5/32.png) [@Stephen](https://discourse.julialang.org/u/Stephen)\
**Post date:** [January 11, 2025, 4:07am UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/15 "2025-01-11T04:07:34Z")

</div>

Is

> [@Reza\_AO](#):
>
> `h_conv = 20`

match with your `T_exact`? As the paper’s Fig. 7. A shows the higher `heq` shifts the curve to left.  
BTW, is x dir represents Radial?

---

<div class="post-metadata">

**Author:** ![Reza\_AO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reza_ao/32/51137_2.png) [@Reza\_AO](https://discourse.julialang.org/u/Reza_AO)\
**Post date:** [January 14, 2025, 12:28pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/16 "2025-01-14T12:28:43Z")

</div>

The problem is not from parameters. I have checked them out. The problem is the delay. I assume the chosen grids and timesteps make my problem stable and consistent (I used small values and there was no change in the results).

> [@Stephen](#):
>
> BTW, is x dir represents Radial?

Yes, as it must be.

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [January 14, 2025, 12:54pm UTC](https://discourse.julialang.org/t/molfinitedifference-for-2d-transient-heat-equation-with-heat-generation/124498/17 "2025-01-14T12:54:49Z")

</div>

If it is radial distance, I think the correct diffusion operator should be

```julia
kₓ*(Dx(x*Dx(T(x,y,t))))/x

```

and you get

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/b/5bfe94ff81a24ba40dec4412e2d555593489b007.png)

But if you used exactly the same equation `kₓ*Dxx(T(x,y,t))` in your reference solution it still does not explain the delay you observed. If on the other hand your reference solution discretized the correct radial operator, then maybe you have an explanation.