# Facing problem writing code of heat equation

**URL:** <https://discourse.julialang.org/t/facing-problem-writing-code-of-heat-equation/127524>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [March 30, 2025, 8:37pm UTC](https://discourse.julialang.org/t/facing-problem-writing-code-of-heat-equation/127524 "2025-03-30T20:37:19Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Stayal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stayal/32/215692_2.png) [@Stayal](https://discourse.julialang.org/u/Stayal)\
**Post date:** [March 30, 2025, 8:37pm UTC](https://discourse.julialang.org/t/facing-problem-writing-code-of-heat-equation/127524/1 "2025-03-30T20:37:19Z")

</div>

Hello everyone

I am a beginner in Gridap. I try to write heat equation code but it shows error

```julia
using Gridap
import Gridap: ∇
using Printf

phi(x) = x[1]*(1 - x[1]) * x[2]*(1 - x[2])
u_ex(x, t) = t^2 * phi(x)
f(x, t) = 2*t*phi(x) + 2*t^2*( x[1]*(1 - x[1]) + x[2]*(1 - x[2]) )

function generate_model_unit_square(nk;simplexify_model=true)
    domain = (0.0,1.0,0.0,1.0)
    n = 2^nk
    #Discrete model
    partition = (n,n)
    model = CartesianDiscreteModel(domain, partition)
    if simplexify_model
      simplexify(model)
    else
      model
    end
end

function run(model; T=1.0, num_steps=256)
  dt = T / num_steps
  order = 1
  reffe = ReferenceFE(lagrangian,Float64,order)
  V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
  U = TrialFESpace(V0,u_ex)

  degree = 2*order
  Ω = Triangulation(model)
  dΩ = Measure(Ω,degree)
  uh_prev = zero(U)
  a(u,v) = (1/dt)*u*v*dΩ + ∫( ∇(v)⊙∇(u) )*dΩ
  for step in 1:num_steps
      t = step * dt
      ft(x) = f(x, t)
      b(v) = ∫( v*ft )*dΩ + (1/dt)*∫( v * uh_prev )*dΩ
      op = AffineFEOperator(a,b,U,V0)
      uh = solve(op)
      uh_prev = uh
      e = u_ex - uh
      el2 = sqrt(sum(∫(e*e)*dΩ))
  el2

end

function conv_test(nkmax)

  error = Float64[]
  hh = Float64[]
  rate = Float64[]
  
  push!(rate,0.)
  for nk in 1:nkmax
      println(" ******** Refinement step: $nk")
      model=generate_model_unit_square(nk;simplexify_model=true)
        push!(hh,sqrt(2)/2^nk)
        el2 =run(model;T=1.0, num_steps=256)

        push!(error,el2)
        if nk>1
            push!(rate, log(error[nk]/error[nk-1])/log(hh[nk]/hh[nk-1]))
        end
    end
    println("========================================================")
    println(" h & e(uS) & r(uS) ")
    println("========================================================")
    for nk in 1:nkmax
        @printf(" %.4f & %.2e & %.3f \n", hh[nk], error[nk], rate[nk], );
    end
    println("========================================================")
end
conv_test(4) 

```

Please help in writing correct code for heat equation.

Thanks  
Shivangi