Facing problem writing code of heat equation

New to Julia
1

Hello everyone

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

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

1 Like