Convergence Error

[Narwal-Pooja]

Hello everyone

When I run the code for the stationary Stokes equation, it works fine; however, when I run the code for the time-dependent Stokes equation, the error remains constant and does not decrease with mesh refinement. I suspect that there might be a syntax error.

using Gridap
import Gridap: ∇
using Printf

u_ex(t) = x -> VectorValue(t*sin(pi*x[1])*sin(pi*x[1])*sin(pi*x[2])*cos(pi*x[2]), -t*sin(pi*x[2])*sin(pi*x[2])*sin(pi*x[1])*cos(pi*x[1]))
p_ex(t) = x -> sin(2*pi*x[1])*sin(2*pi*x[2])
uu_ex(t) = x -> VectorValue(sin(pi*x[1])*sin(pi*x[1])*sin(pi*x[2])*cos(pi*x[2]), -sin(pi*x[2])*sin(pi*x[2])*sin(pi*x[1])*cos(pi*x[1]))

function 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)
  order = 2 
  reffeᵤ  = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
  reffeₚ = ReferenceFE(lagrangian,Float64,order-1)
  V0 = TestFESpace(model, reffeᵤ, conformity=:H1, dirichlet_tags=["boundary"])
  Q = TestFESpace(model, reffeₚ, conformity=:C0, constraint=:zeromean)
  U = TransientTrialFESpace(V0,[u_ex])
  P = TrialFESpace(Q)
  
  X = MultiFieldFESpace([U,P])
  Y = MultiFieldFESpace([V0,Q])
  
  degree = 2*order
  Ω = Triangulation(model)
  dΩ = Measure(Ω,degree)

  f(t) = x -> uu_ex(t)(x) - (∇⋅∇(u_ex(t)))(x) + (∇(p_ex(t)))(x)
  g(t) = x -> (∇⋅u_ex(t))(x)

  m(t,(u,p), (v,q)) = ∫(v ⋅ u)dΩ
  L(t,(v,q))= ∫(v⋅f(t))dΩ + ∫(q*g(t))dΩ
  A(t,(u,p),(v,q)) = ∫(∇(u)⊙∇(v))dΩ + ∫(q*(∇⋅u))dΩ - ∫(p*(∇⋅v))dΩ 

  op = TransientLinearFEOperator((A,m), L, X, Y)
  ls = LUSolver()
  Δt = 0.01
  θ = 1
  solver = ThetaMethod(ls, Δt, θ)
  tableau = :SDIRK_2_2
  solver_rk = RungeKutta(ls, ls, Δt, tableau)
  t0, tF = 0.0, 1.0

  xh0 = interpolate_everywhere([u_ex(t0),p_ex(t0)],X(t0))
  uh, ph = solve(solver, op, t0, tF, xh0)
  
  uh_array = collect(uh)
 
  uh_final = uh_array[end][2] 

  exact_final = u_ex(tF)
  
  e = uh_final - exact_final
  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=model_unit_square(nk;simplexify_model=true)
        push!(hh,sqrt(2)/2^nk)
        el2 =run(model)

        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(5)

RESULT

******** Refinement step: 1
******** Refinement step: 2
******** Refinement step: 3
******** Refinement step: 4
******** Refinement step: 5
===================================
     h   &   e(uS)  & r(uS)  
===================================
 0.7071 & 7.07e-01 & 0.000  
 0.3536 & 7.63e-01 & -0.109  
 0.1768 & 7.70e-01 & -0.014  
 0.0884 & 7.71e-01 & -0.001  
 0.0442 & 7.71e-01 & -0.000  
===================================

[Narwal-Pooja]

Please provide some suggestions

[JordiManyer]

I am quite sure you have not found this in any tutorial:

  uh, ph = solve(solver, op, t0, tF, xh0)
  uh_array = collect(uh)
  uh_final = uh_array[end][2] 

Although solve returns an iterable, caches are reused throughout. This means there is aliasing between different iterations, otherwise the memory costs would be prohibitive. Also the solve method returns (t,x), not (u,p). Unless you go into the source code and develop an understanding of how things work, I would try to follow the tutorials/tests.

fesltn = solve(odeslvr, tfeop, t0, tF, xhs0)
for (t_n, xhs_n) in fesltn
  eh_u = u(t_n) - xhs_n[1]
  e_u = sqrt(sum(∫(eh_n ⋅ eh_n) * dΩ))

  eh_p = p(t_n) - xhs_n[2]
  e_p = sqrt(sum(∫(eh_n ⋅ eh_n) * dΩ))
end

[Narwal-Pooja]

Thank you