In this tutorial we will learn how to use
GridapODEs.jl for approximating transient PDEs by using time marching schemes (method of lines). We consider the heat equation, a.k.a. the transient Poisson equation.
We will focus on the time discretization on the equations, assuming that the reader is familiar with the
Gridap API for spatial finite element discretizations. See, e.g., tutorial 1 for more details.
We solve the heat equation in a 2-dimensional domain $\Omega$, the unit square, with Homogeneous Dirichlet boundaries on the whole boundary $\partial \Omega$. We consider a time-dependent conductivity $\kappa(t)=1.0 + 0.95\sin(2\pi t)$, a time-dependent volumetric forcing term $f(t) = \sin(\pi t)$ and a constant Homogeneous boundary condition $g = 0.0$. The initial solution is $u(x,0) = u_0 = 0$. With these definitions, the strong form of the problem reads:
The weak form of the problem reads: find $u(t)\in U_g(t)$ such that
Note that $U_g(t)$ is a transient FE space, in the sense that Dirichlet boundary value of functions in $U_g$ can change in time (even though this is not the case in this tutorial). The definition of $m(u,v)$, $a(u,v)$ and $b(v)$ is as follows:
As usual, let us first load
First, we define the
DiscreteModel and the
Triangulation. More details on this can be found in tutorial 2.
𝒯 = CartesianDiscreteModel((0,1,0,1),(20,20)) Ω = Interior(𝒯) dΩ = Measure(Ω,2)
In this tutorial we will use linear Lagrangian Finite Elements.
refFE = ReferenceFE(lagrangian,Float64,1)
The space of test functions is constant in time and is defined in steady problems:
V = TestFESpace(𝒯,refFE,dirichlet_tags="boundary")
The trial space is now a
TransientTrialFESpace, which is constructed from a
TestFESpace and a function (or vector of functions) for the Dirichlet boundary condition/s. In that case, the boundary condition function is a time-independent constant, but it could also be a time-dependent field depending on the coordinates $x$ and time $t$.
g(x,t::Real) = 0.0 g(t::Real) = x -> g(x,t) U = TransientTrialFESpace(V,g)
The weak form of the problem follows the same structure as other
Gridap tutorials, where we define the bilinear and linear forms to define the
FEOperator. In this case we need to deal with time-dependent quantities and with the presence of time derivatives. The former is handled by passing the time, $t$, as an additional argument to the form, i.e. $a(t,u,v)$. The latter is defined using the time derivative operator
The most general way of constructing a transient FE operator is by using the
TransientFEOperator function, which receives a residual, a Jacobian with respect to the unknown and a Jacobian with respect to the time derivative.
κ(t) = 1.0 + 0.95*sin(2π*t) f(t) = sin(π*t) res(t,u,v) = ∫( ∂t(u)*v + κ(t)*(∇(u)⋅∇(v)) - f(t)*v )dΩ jac(t,u,du,v) = ∫( κ(t)*(∇(du)⋅∇(v)) )dΩ jac_t(t,u,duₜ,v) = ∫( duₜ*v )dΩ op = TransientFEOperator(res,jac,jac_t,U,V)
We can also take advantage of automatic differentiation techniques to compute both Jacobians and use the
TransientFEOperator function sending just the residual.
op_AD = TransientFEOperator(res,U,V)
Alternatively, we can exploit the fact that the problem is linear and use the transient Affine FE operator signature
TransientAffineFEOperator. In that case, we send a form for the mass contribution, $m$, a form for the stiffness contribution, $a$, and the forcing term, $b$.
m(t,u,v) = ∫( u*v )dΩ a(t,u,v) = ∫( κ(t)*(∇(u)⋅∇(v)) )dΩ b(t,v) = ∫( f(t)*v )dΩ op_Af = TransientAffineFEOperator(m,a,b,U,V)
For time-dependent problems with constant coefficients, which is not the case of this tutorial, one could use the optimized operator
TransientConstantMatrixFEOperator, which assumes that the matrix contributions ($m$ and $a$) are time-independent. That is:
m₀(u,v) = ∫( u*v )dΩ a₀(u,v) = ∫( κ(0.0)*(∇(u)⋅∇(v)) )dΩ op_CM = TransientConstantMatrixFEOperator(m₀,a₀,b,U,V)
Going further, if we had a problem with constant forcing term, i.e. constant force and constant boundary conditions, we could have used the
TransientConstantFEOperator. In that case the linear form is also time-independent.
b₀(v) = ∫( f(0.0)*v )dΩ op_C = TransientConstantFEOperator(m₀,a₀,b₀,U,V)
Once we have the FE operator defined, we proceed with the definition of the transient solver. First, we define a linear solver to be used at each time step. Here we use the
LUSolver, but other choices are possible.
linear_solver = LUSolver()
Then, we define the ODE solver. That is, the scheme that will be used for the time integration. In this tutorial we use the
ThetaMethod with $\theta = 0.5$, resulting in a 2nd order scheme. The
ThetaMethod function receives the linear solver, the time step size $\Delta t$ (constant) and the value of $\theta$.
Δt = 0.05 θ = 0.5 ode_solver = ThetaMethod(linear_solver,Δt,θ)
Finally, we define the solution using the
solve function, giving the ODE solver, the FE operator, an initial solution, an initial time and a final time. To construct the initial condition we interpolate the initial value (in that case a constant value of 0.0) into the FE space $U(t)$ at $t=0.0$.
u₀ = interpolate_everywhere(0.0,U(0.0)) t₀ = 0.0 T = 10.0 uₕₜ = solve(ode_solver,op,u₀,t₀,T)
We should highlight that
uₕₜ is just an iterable function and the results at each time steps are only computed when iterating over it, i.e., lazily. We can post-process the results and generate the corresponding
vtk files using the
createvtk functions. The former will create a
.pvd file with the collection of
.vtu files saved at each time step by
createvtk. The computation of the problem solutions will be triggered in the following loop:
createpvd("poisson_transient_solution") do pvd for (uₕ,t) in uₕₜ pvd[t] = createvtk(Ω,"poisson_transient_solution_$t"*".vtu",cellfields=["u"=>uₕ]) end end
This page was generated using Literate.jl.