Hi all-
Is it currently possible to solve transient problems for multifield FE spaces using Gridap.jl? I’m currently trying to solve a Boussinesq Navier-Stokes system, and am running into some issues, so I’d like to know if it’s at least possible with the native current release before getting too much more invested.
If not, I’ll turn to Trixi.jl, where I believe this is possible.
My current code is below. I haven’t put any thought into the ODE solver yet: just trying to get this off the ground. My interest in using Gridap.jl/Trixi.jl is premised on the hope that it will be easy to modify the physics and iterate on the system being solved.
using Gridap, GridapGmsh, Gridap.Io
# inputs: mesh file, initial condition, parameters for forcing function
# parameters
ic = [VectorValue(0, 0, 0), 0, 80.0] # read in the initial condition from an input file
day = 1 # read in the day from an input file
h₀ = 70.0 # mean forcing on Dirichlet boundary
h₁ = 10.0 # amplitude of the forcing on Dirichlet boundary
const ℛ = 10.0 # Reynolds number
const 𝒫 = 10.0 # Peclet number
# We load in the MESH..
msh_file = joinpath(dirname(@__DIR__),"system/box-1e4.msh")
model = GmshDiscreteModel(msh_file)
# .. and indicate the DIRICHLET BOUNDARIES:
dirichlet_tags = [
"westFace",
"eastFace",
"southFace",
"northFace",
# "westBase",
# "eastBase",
# "southBase",
# "northBase",
"upFace",
"bottom"
]
dirichlet_tags_velocity = dirichlet_tags
dirichlet_tags_temperature = dirichlet_tags
# .. then, selecting the reference order for the lagrange elements, ..
order = 2
# .. we define the REFERENCE FINITE ELEMENTS:
reffe_u = ReferenceFE(lagrangian,VectorValue{3,Float64},order) # .... of velocity, reffe_u,
reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) # .... of pressure, reffe_p,
reffe_θ = ReferenceFE(lagrangian,Float64,order) # .... and of temperature, reffe_θ.
# ... from which we build the TEST SPACES:
V = TestFESpace(model,reffe_u,conformity=:H1,dirichlet_tags=dirichlet_tags_velocity) # .... of velocity, V,
Q = TestFESpace(model,reffe_p,conformity=:L2,constraint=:zeromean) # .... of pressure, Q,
Φ = TestFESpace(model,reffe_θ,conformity=:H1,dirichlet_tags=dirichlet_tags_temperature) # .... and of temperature, Φ.
degree = 2 # degree of the quadrature rule
Ω = Triangulation(model) # for CartesianDiscreteModel, replace with Interior(model)
dΩ = Measure(Ω,degree) # approximation of Lebesgue measure
Γ = BoundaryTriangulation(model) # for CartesianDiscreteModel, replace with Boundary(model)
dΓ = Measure(Γ,degree)
# transient problem: https://gridap.github.io/Tutorials/dev/pages/t017_poisson_transient/
# (steady) Navier-Stokes equations: https://gridap.github.io/Tutorials/dev/pages/t008_inc_navier_stokes/
h(x,t::Real) = h₁ * sin(t/24.0 * 2*π) + h₀ # point to a function that returns the value of the forcing function
h(t::Real) = x -> h(x,t)
# ... from which we build the TRIAL SPACES:
uΓ = VectorValue(0.0,0.0,0.0) # forcing for velocity
U = TransientTrialFESpace(V,x->uΓ) # .... of velocity, U,
P = TransientTrialFESpace(Q,x->0.0) # .... of pressure, P,
Θh = TransientTrialFESpace(Φ,h) # .... and of temperature, Θh, which is transient, depending on the Dirichlet data, h(x,t).
# ... which we wrap us as a MULTI-FIELD SPACE:
X = TransientMultiFieldFESpace([U,P,Θh])
Y = MultiFieldFESpace([V,Q,Φ])
# [dev] runs up to here
# To build the weak form, we first define a convection operator
conv(u,∇u) = ℛ*(∇u')⋅u
dconv(du,∇du,u,∇u) = conv(u,∇du)+conv(du,∇u)
# Define j, the unit vector in the z-direction
j = [0.0,0.0,1.0]
# Define the weak form of the transient problem through specificication of...
m((u,θ),(v,ϕ)) = ∫( ∂t(u)⋅v )dΩ + ∫( ∂t(θ)*ϕ )dΩ # ... a transient term,
a((u,p,θ),(v,q,ϕ)) = ν*∫( ∇(u)⊙∇(v) )dΩ + ϵ*∫( ∇(θ)⋅∇(ϕ) )dΩ + ∫( q*(∇⋅u) )dΩ - ∫( p*(∇⋅v) )dΩ - ∫( θ*(j⋅v) )dΩ # ... the linear stiffness components (⊙ is Hadamard product),
c((u,θ),(v,ϕ)) = ∫( v⊙(conv∘(u,∇(u))) )dΩ + ∫( ϕ*(u⋅∇(θ)) )dΩ# ... the convective (nonlinear) components (∘ is composition),
b(t,v) = 0.0 # ... and the forcing vector.
# We combine these terms into the following residual:
res((u,p,θ),(v,q,ϕ)) = m((u,θ),(v,ϕ)) + a((u,p,θ),(v,q,ϕ)) + c(u,v) - b(t,v)
# To define the jacobian, we must specify the derivative of the nonlinear terms:
dc((u,θ),(du,dθ),(v,ϕ)) = ∫( v⊙(dconv∘(du,∇(du),u,∇(u))) )dΩ + ∫( ϕ*(du⋅∇(θ)+u⋅∇(dθ)) )dΩ
# and thus our jacobian is given by
jac((u,p,θ),(du,dp,dθ),(v,q,ϕ)) = a((du,dp,dθ),(v,q,ϕ)) + dc((u,θ),(du,dθ),(v,ϕ))
# and the time derivative of the jacobian is given by
jac_t((u,p,θ),(du,dp,dθ),(v,q,ϕ)) = m((du,dθ),(v,ϕ))
op = TransientFEOperator(res,jac,jac_t,X,Y)
# generate outputs at multiple time steps
linear_solver = LUSolver()
Δt = 0.1
ode_solver = BackwardEuler(linear_solver,Δt)
# interpolate the initial condition
x₀ = interpolate_everywhere(ic,X(0.0)) # needs to match the DOF of the space
t₀ = 0.0
T = 1.0
# https://www.youtube.com/watch?v=heeiSoKnlUk
xₕₜ = solve(ode_solver,op,x₀,t₀,T)
# include pointer for where to save the outputs
# this isn't working, need to be patient and figure out where to save the outputs
createpvd("building-twins/model-1/VTK/boussinesq_solution") do pvd
for (xₕ,t) in xₕₜ
pvd[t] = createvtk(Ω,"building-twins/model-1/VTK/boussinesq_solution_$(round(t,digits=3))"*".vtu",cellfields=["u"=>xₕ[1],"p"=>xₕ[2],"θ"=>xₕ[3]])
end
end