Facing problem in time dependent Stokes code

[Narwal-Pooja]

Hello

Below I attached both the error and code

CODE:

using Gridap
import Gridap: ∇
using Printf


u_ex(t) = x -> VectorValue(x[1],x[2])*t
p_ex(t) = x ->  (x[1]-x[2])*t
uu_ex(t) = x -> VectorValue(x[1],x[2])
domain = (0.0,1.0,0.0,1.0)
n      = 32
partition = (n,n)
model = CartesianDiscreteModel(domain, partition)
order = 2
reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
V0 = TestFESpace(model, reffeᵤ, conformity=:H1, dirichlet_tags="boundary")
reffeₚ = ReferenceFE(lagrangian,Float64,order-1)
Q = TestFESpace(model, reffeₚ, conformity=:L2, constraint=:zeromean)

U = TransientTrialFESpace(V0,u_ex)
P = TrialFESpace(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)

X = MultiFieldFESpace([U,P])
Y = MultiFieldFESpace([V0,Q])

m(t,dtu,v) = ∫(v * dtu)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

uh0 = interpolate_everywhere(u_ex(t0),U(t0))
ph0 = interpolate_everywhere(p_ex(t0),P(t0))
xh0 = interpolate_everywhere([uh0,ph0],X(t0))

uh, ph = solve(solver, op, t0, tF, xh0)

ERROR :

pooja@pooja-Inspiron-3501:~$ julia --project=~/gridap_project ~/Downloads/julia/stokes.jl
ERROR: LoadError: Function get_data is not implemented for MultiFieldCellField at this moment.
You need to extract the individual fields and then evaluate them separately.

If ever implement this, evaluating a `MultiFieldCellField` directly would provide,
at each evaluation point, a tuple with the value of the different fields.

Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] macro expansion
    @ ~/.julia/packages/Gridap/lHC3y/src/Helpers/Macros.jl:21 [inlined]
  [3] get_data(f::Gridap.MultiField.MultiFieldCellField{ReferenceDomain})
    @ Gridap.MultiField ~/.julia/packages/Gridap/lHC3y/src/MultiField/MultiFieldCellFields.jl:28
  [4] evaluate!(cache::Nothing, f::Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, x::Gridap.CellData.CellPoint{ReferenceDomain, FillArrays.Fill{Vector{VectorValue{2, Float64}}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Gridap.Geometry.CartesianCoordinates{2, Float64, typeof(identity)}}}, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Vector{VectorValue{2, Float64}}, 2, Tuple{Gridap.Geometry.CartesianCellNodes{2}}}, Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}, Gridap.Geometry.Injective}})
    @ Gridap.CellData ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:401
  [5] evaluate
    @ ~/.julia/packages/Gridap/lHC3y/src/Arrays/Maps.jl:87 [inlined]
  [6] CellField
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:397 [inlined]
  [7] #11
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:495 [inlined]
  [8] map(f::Gridap.CellData.var"#11#13"{Gridap.CellData.CellPoint{ReferenceDomain, FillArrays.Fill{Vector{VectorValue{2, Float64}}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Gridap.Geometry.CartesianCoordinates{2, Float64, typeof(identity)}}}, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Vector{VectorValue{2, Float64}}, 2, Tuple{Gridap.Geometry.CartesianCellNodes{2}}}, Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}, Gridap.Geometry.Injective}}}, t::Tuple{Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, Gridap.MultiField.MultiFieldCellField{ReferenceDomain}})
    @ Base ./tuple.jl:356
  [9] macro expansion
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:495 [inlined]
 [10] macro expansion
    @ ~/.julia/packages/Gridap/lHC3y/src/Helpers/Macros.jl:61 [inlined]
 [11] Gridap.CellData.OperationCellField(::Operation{typeof(*)}, ::Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, ::Vararg{Gridap.MultiField.MultiFieldCellField{ReferenceDomain}})
    @ Gridap.CellData ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:489
 [12] _operate_cellfields(::Operation{typeof(*)}, ::Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, ::Vararg{Gridap.MultiField.MultiFieldCellField{ReferenceDomain}})
    @ Gridap.CellData ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:558
 [13] evaluate!
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:451 [inlined]
 [14] evaluate
    @ ~/.julia/packages/Gridap/lHC3y/src/Arrays/Maps.jl:87 [inlined]
 [15] Map
    @ ~/.julia/packages/Gridap/lHC3y/src/Arrays/Maps.jl:90 [inlined]
 [16] *
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:638 [inlined]
 [17] m(t::Float64, dtu::Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, v::Gridap.MultiField.MultiFieldCellField{ReferenceDomain})
    @ Main ~/Downloads/julia/stokes.jl:32
 [18] (::Gridap.ODEs.var"#112#116"{Int64, Tuple{typeof(A), typeof(m)}})(t::Float64, u::Gridap.ODEs.TransientMultiFieldCellField{Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{ReferenceDomain}}}, duk::Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, v::Gridap.MultiField.MultiFieldCellField{ReferenceDomain})
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/TransientFEOperators.jl:693
 [19] allocate_odeopcache(odeop::Gridap.ODEs.ODEOpFromTFEOp{Gridap.ODEs.LinearODE}, t::Float64, us::Tuple{Vector{Float64}, Vector{Float64}})
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/ODEOpsFromTFEOps.jl:109
 [20] allocate_odecache(odeslvr::ThetaMethod, odeop::Gridap.ODEs.ODEOpFromTFEOp{Gridap.ODEs.LinearODE}, t0::Float64, us0::Tuple{Vector{Float64}})
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/ODESolvers/ThetaMethod.jl:116
 [21] iterate(odesltn::Gridap.ODEs.GenericODESolution)
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/ODESolutions.jl:60
 [22] iterate(tfesltn::Gridap.ODEs.GenericTransientFESolution)
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/TransientFESolutions.jl:75
 [23] indexed_iterate(I::Gridap.ODEs.GenericTransientFESolution, i::Int64)
    @ Base ./tuple.jl:162
 [24] top-level scope
    @ ~/Downloads/julia/stokes.jl:49
in expression starting at /home/pooja/Downloads/julia/stokes.jl:49

[JordiManyer]

The issue is in m, where v should be (v,q). I.e

m(t,dtu,(v,q)) = ∫(v * dtu)dΩ

[Narwal-Pooja]

Thank you for the response but I am still getting error attached below :

julia --project=~/gridap_project ~/Downloads/julia/stokes.jl
ERROR: LoadError: Function get_data is not implemented for MultiFieldCellField at this moment.
You need to extract the individual fields and then evaluate them separately.

If ever implement this, evaluating a `MultiFieldCellField` directly would provide,
at each evaluation point, a tuple with the value of the different fields.

Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] macro expansion
    @ ~/.julia/packages/Gridap/lHC3y/src/Helpers/Macros.jl:21 [inlined]
  [3] get_data(f::Gridap.MultiField.MultiFieldCellField{ReferenceDomain})
    @ Gridap.MultiField ~/.julia/packages/Gridap/lHC3y/src/MultiField/MultiFieldCellFields.jl:28
  [4] evaluate!(cache::Nothing, f::Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, x::Gridap.CellData.CellPoint{ReferenceDomain, FillArrays.Fill{Vector{VectorValue{2, Float64}}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Gridap.Geometry.CartesianCoordinates{2, Float64, typeof(identity)}}}, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Vector{VectorValue{2, Float64}}, 2, Tuple{Gridap.Geometry.CartesianCellNodes{2}}}, Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}, Gridap.Geometry.Injective}})
    @ Gridap.CellData ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:401
  [5] evaluate
    @ ~/.julia/packages/Gridap/lHC3y/src/Arrays/Maps.jl:87 [inlined]
  [6] CellField
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:397 [inlined]
  [7] #11
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:495 [inlined]
  [8] map(f::Gridap.CellData.var"#11#13"{Gridap.CellData.CellPoint{ReferenceDomain, FillArrays.Fill{Vector{VectorValue{2, Float64}}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Gridap.Geometry.CartesianCoordinates{2, Float64, typeof(identity)}}}, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Vector{VectorValue{2, Float64}}, 2, Tuple{Gridap.Geometry.CartesianCellNodes{2}}}, Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}, Gridap.Geometry.Injective}}}, t::Tuple{Gridap.MultiField.MultiFieldFEBasisComponent{Gridap.FESpaces.SingleFieldFEBasis{Gridap.FESpaces.TestBasis, ReferenceDomain}}, Gridap.MultiField.MultiFieldCellField{ReferenceDomain}})
    @ Base ./tuple.jl:356
  [9] macro expansion
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:495 [inlined]
 [10] macro expansion
    @ ~/.julia/packages/Gridap/lHC3y/src/Helpers/Macros.jl:61 [inlined]
 [11] Gridap.CellData.OperationCellField(::Operation{typeof(*)}, ::Gridap.MultiField.MultiFieldFEBasisComponent{Gridap.FESpaces.SingleFieldFEBasis{Gridap.FESpaces.TestBasis, ReferenceDomain}}, ::Vararg{CellField})
    @ Gridap.CellData ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:489
 [12] _operate_cellfields(::Operation{typeof(*)}, ::Gridap.MultiField.MultiFieldFEBasisComponent{Gridap.FESpaces.SingleFieldFEBasis{Gridap.FESpaces.TestBasis, ReferenceDomain}}, ::Vararg{Any})
    @ Gridap.CellData ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:558
 [13] evaluate!
    @ ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:451 [inlined]
 [14] evaluate
    @ ~/.julia/packages/Gridap/lHC3y/src/Arrays/Maps.jl:87 [inlined]
 [15] Map
    @ ~/.julia/packages/Gridap/lHC3y/src/Arrays/Maps.jl:90 [inlined]
 [16] *(a::Gridap.MultiField.MultiFieldFEBasisComponent{Gridap.FESpaces.SingleFieldFEBasis{Gridap.FESpaces.TestBasis, ReferenceDomain}}, b::Gridap.MultiField.MultiFieldCellField{ReferenceDomain})
    @ Gridap.CellData ~/.julia/packages/Gridap/lHC3y/src/CellData/CellFields.jl:638
 [17] m(t::Float64, dtu::Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, ::Gridap.MultiField.MultiFieldCellField{ReferenceDomain})
    @ Main ~/Downloads/julia/stokes.jl:32
 [18] (::Gridap.ODEs.var"#112#116"{Int64, Tuple{typeof(A), typeof(m)}})(t::Float64, u::Gridap.ODEs.TransientMultiFieldCellField{Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{ReferenceDomain}}}, duk::Gridap.MultiField.MultiFieldCellField{ReferenceDomain}, v::Gridap.MultiField.MultiFieldCellField{ReferenceDomain})
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/TransientFEOperators.jl:693
 [19] allocate_odeopcache(odeop::Gridap.ODEs.ODEOpFromTFEOp{Gridap.ODEs.LinearODE}, t::Float64, us::Tuple{Vector{Float64}, Vector{Float64}})
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/ODEOpsFromTFEOps.jl:109
 [20] allocate_odecache(odeslvr::ThetaMethod, odeop::Gridap.ODEs.ODEOpFromTFEOp{Gridap.ODEs.LinearODE}, t0::Float64, us0::Tuple{Vector{Float64}})
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/ODESolvers/ThetaMethod.jl:116
 [21] iterate(odesltn::Gridap.ODEs.GenericODESolution)
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/ODESolutions.jl:60
 [22] iterate(tfesltn::Gridap.ODEs.GenericTransientFESolution)
    @ Gridap.ODEs ~/.julia/packages/Gridap/lHC3y/src/ODEs/TransientFESolutions.jl:75
 [23] indexed_iterate(I::Gridap.ODEs.GenericTransientFESolution, i::Int64)
    @ Base ./tuple.jl:162
 [24] top-level scope
    @ ~/Downloads/julia/stokes.jl:49
in expression starting at /home/pooja/Downloads/julia/stokes.jl:49
pooja@pooja-Inspiron-3501:~$ 

[JordiManyer]

Same for the trial variable.

m(t,(u,p),(v,q)) = ∫(v * u)dΩ

the error is telling you that you are using a MultiFieldCellField in a non-intended way. This is because your trial/test variables contain more than one variable, but you are treating them as one.

[Narwal-Pooja]

Thank you. I tried with both m(t,(u,p),(v,q)) = ∫(v * dtu)dΩ or m(t,(dtu,p),(v,q)) = ∫(v * dtu)dΩ but doesn't work.

[JordiManyer]

Can you post the code and error when changing to

m(t,(u,p),(v,q)) = ∫(v * u)dΩ

[Narwal-Pooja]

Sir, if we define

  m(t, (u, p), (v, q)) = ∫ (v · u) dΩ,

how does the time term work in the code? In the transient Poisson code, we use

  m(t, dtu, v) = ∫ (v · dtu) dΩ to handle du/dt term, That's why I use m(t,(u,p),(v,q)) = ∫(v * dtu)dΩ in Stokes code.

Code

using Gridap
import Gridap: ∇
using Printf


u_ex(t) = x -> VectorValue(x[1],x[2])*t
p_ex(t) = x ->  (x[1]-x[2])*t
uu_ex(t) = x -> VectorValue(x[1],x[2])
domain = (0.0,1.0,0.0,1.0)
n      = 32
partition = (n,n)
model = CartesianDiscreteModel(domain, partition)
order = 2
reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
V0 = TestFESpace(model, reffeᵤ, conformity=:H1, dirichlet_tags="boundary")
reffeₚ = ReferenceFE(lagrangian,Float64,order-1)
Q = TestFESpace(model, reffeₚ, conformity=:L2, constraint=:zeromean)

U = TransientTrialFESpace(V0,u_ex)
P = TrialFESpace(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)

X = MultiFieldFESpace([U,P])
Y = MultiFieldFESpace([V0,Q])

m(t,(u,p),(v,q)) = ∫(v * dtu)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

uh0 = interpolate_everywhere(u_ex(t0),U(t0))
ph0 = interpolate_everywhere(p_ex(t0),P(t0))
xh0 = interpolate_everywhere([uh0,ph0],X(t0))

uh, ph = solve(solver, op, t0, tF, xh0)

[JordiManyer]

They are just variables, the name doesnt matter. In your function, dtu is not defined.

As far as I understand, in transient stokes there is only a time derivative on the velocity, not the pressure. So you only need a mass matrix on (u,v).

[JordiManyer]

https://github.com/gridap/Gridap.jl/blob/master/test%2FODEsTests%2FTransientFEProblemsTests%2FStokesEquationTests.jl

[Narwal-Pooja]

Thank you.