Merge pull request #63 from gridap/new_API_transient_fe_operators

New api transient fe operators

README.md

@@ -14,9 +14,6 @@ using ForwardDiff
 using GridapODEs.ODETools
 using GridapODEs.TransientFETools
 
-import Gridap: ∇
-import GridapODEs.TransientFETools: ∂t
-
 θ = 0.5
 
 u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
@@ -48,9 +45,9 @@ a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ
 b(v,t) = ∫( v*f(t) )dΩ
 m(u,v) = ∫( v*u )dΩ
 
-res(t,(u,ut),v) = a(u,v) + m(ut,v) - b(v,t)
-jac(t,(u,ut),du,v) = a(du,v)
-jac_t(t,(u,ut),dut,v) = m(dut,v)
+res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t)
+jac(t,u,du,v) = a(du,v)
+jac_t(t,u,dut,v) = m(dut,v)
 
 op = TransientFEOperator(res,jac,jac_t,U,V0)
 

src/TransientFETools/ODEOperatorInterfaces.jl

@@ -61,10 +61,11 @@ function residual!(
   xhF::Tuple{Vararg{AbstractVector}},
   ode_cache)
   Xh, = ode_cache
-  xh = ()
-  for i in 1:get_order(op)+1
-    xh = (xh...,EvaluationFunction(Xh[i],xhF[i]))
+  dxh = ()
+  for i in 2:get_order(op)+1
+    dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i]))
   end
+  xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh)
   residual!(b,op.feop,t,xh,ode_cache)
 end
 
@@ -86,10 +87,11 @@ function jacobian!(
   γᵢ::Real,
   ode_cache)
   Xh, = ode_cache
-  xh = ()
-  for i in 1:get_order(op)+1
-    xh = (xh...,EvaluationFunction(Xh[i],xhF[i]))
+  dxh = ()
+  for i in 2:get_order(op)+1
+    dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i]))
   end
+  xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh)
   jacobian!(A,op.feop,t,xh,i,γᵢ,ode_cache)
 end
 
@@ -104,9 +106,10 @@ function jacobians!(
   γ::Tuple{Vararg{Real}},
   ode_cache)
   Xh, = ode_cache
-  xh = ()
-  for i in 1:get_order(op)+1
-    xh = (xh...,EvaluationFunction(Xh[i],xhF[i]))
+  dxh = ()
+  for i in 2:get_order(op)+1
+    dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i]))
   end
+  xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh)
   jacobians!(J,op.feop,t,xh,γ,ode_cache)
 end

src/TransientFETools/TransientCellField.jl

@@ -0,0 +1,69 @@
+# Transient CellField
+struct TransientCellField{A} <: CellField
+  cellfield::A
+  derivatives::Tuple#{Vararg{A,B} where B}
+end
+
+# CellField methods
+get_data(f::TransientCellField) = get_data(f.cellfield)
+get_triangulation(f::TransientCellField) = get_triangulation(f.cellfield)
+DomainStyle(::Type{<:TransientCellField{A}}) where A = DomainStyle(A)
+gradient(f::TransientCellField) = gradient(f.cellfield)
+∇∇(f::TransientCellField) = ∇∇(f.cellfield)
+change_domain(f::TransientCellField,trian::Triangulation,target_domain::DomainStyle) = change_domain(f.cellfield,trian,target_domain)
+
+# MultiFieldCellField methods
+num_fields(f::TransientCellField{A}) where A = num_fields(f.cellfield)
+function Base.getindex(f::TransientCellField{A},i::Integer) where A 
+  singleCellfield = getindex(f.cellfield,i)
+  singleDerivatives = ()
+  for i_derivatives in f.derivatives
+    singleDerivatives = (singleDerivatives...,getindex(i_derivatives,i))
+  end
+  # singleDerivatives = (getindex(i_derivatives,i) for i_derivatives in f.derivatives)
+  TransientCellField(singleCellfield,singleDerivatives)
+end
+function Base.iterate(f::TransientCellField{A}) where A
+  cellField, state1 = iterate(f.cellfield)
+  derivatives = ()
+  state2 = []
+  for i_derivatives in f.derivatives
+    single_derivative, single_state = iterate(i_derivatives)
+    derivatives = (derivatives...,single_derivative)
+    push!(state2,single_state)
+  end
+  TransientCellField(cellField,derivatives),(state1,state2)
+end
+function Base.iterate(f::TransientCellField{A},state) where A
+  state1, state2 = state
+  cellField, state1 = iterate(f.cellfield,state[1])
+  derivatives = ()
+  for (i,i_derivatives) in enumerate(f.derivatives)
+    single_derivative, state2[i] = iterate(i_derivatives)
+    derivatives = (derivatives...,single_derivative)
+  end
+  TransientCellField(cellField,derivatives),(state1, state2)
+end
+
+# Transient FEBasis
+struct TransientFEBasis{A} <: FEBasis
+  febasis::A
+  derivatives::Tuple{Vararg{A}}
+end
+
+# FEBasis methods
+get_data(f::TransientFEBasis) = get_data(f.febasis)
+get_triangulation(f::TransientFEBasis) = get_triangulation(f.febasis)
+DomainStyle(::Type{<:TransientFEBasis{A}}) where A = DomainStyle(A)
+BasisStyle(::Type{<:TransientFEBasis{A}}) where A = BasisStyle(A)
+gradient(f::TransientFEBasis) = gradient(f.febasis)
+∇∇(f::TransientFEBasis) = ∇∇(f.febasis)
+change_domain(f::TransientFEBasis,trian::Triangulation,target_domain::DomainStyle) = change_domain(f.febasis,trian,target_domain)
+
+# Time derivative
+function ∂t(f::Union{TransientCellField,TransientFEBasis})
+  cellfield, derivatives = first_and_tail(f.derivatives)
+  TransientCellField(cellfield,derivatives)
+end
+
+∂tt(f::Union{TransientCellField,TransientFEBasis}) = ∂t(∂t(f::Union{TransientCellField,TransientFEBasis}))

src/TransientFETools/TransientFEOperators.jl

@@ -112,27 +112,27 @@ end
 
 function TransientConstantFEOperator(m::Function,a::Function,b::Function,
   trial,test)
-  res(t,(u,ut),v) = -1.0 * b(v)
-  jac(t,(u,ut),du,v) = a(du,v)
-  jac_t(t,(u,ut),dut,v) = m(dut,v)
+  res(t,u,v) = -1.0 * b(v)
+  jac(t,u,du,v) = a(du,v)
+  jac_t(t,u,dut,v) = m(dut,v)
   assem_t = SparseMatrixAssembler(trial,test)
   TransientFEOperatorFromWeakForm{Constant}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1)
 end
 
 function TransientConstantMatrixFEOperator(m::Function,a::Function,b::Function,
   trial,test)
-  res(t,(u,ut),v) = m(ut,v) + a(u,v) - b(t,v)
-  jac(t,(u,ut),du,v) = a(du,v)
-  jac_t(t,(u,ut),dut,v) = m(dut,v)
+  res(t,u,v) = m(∂t(u),v) + a(u,v) - b(t,v)
+  jac(t,u,du,v) = a(du,v)
+  jac_t(t,u,dut,v) = m(dut,v)
   assem_t = SparseMatrixAssembler(trial,test)
   TransientFEOperatorFromWeakForm{ConstantMatrix}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1)
 end
 
 function TransientAffineFEOperator(m::Function,a::Function,b::Function,
   trial,test)
-  res(t,(u,ut),v) = m(t,ut,v) + a(t,u,v) - b(t,v)
-  jac(t,(u,ut),du,v) = a(t,du,v)
-  jac_t(t,(u,ut),dut,v) = m(t,dut,v)
+  res(t,u,v) = m(t,∂t(u),v) + a(t,u,v) - b(t,v)
+  jac(t,u,du,v) = a(t,du,v)
+  jac_t(t,u,dut,v) = m(t,dut,v)
   assem_t = SparseMatrixAssembler(trial,test)
   TransientFEOperatorFromWeakForm{Affine}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1)
 end
@@ -146,10 +146,10 @@ end
 
 function TransientConstantFEOperator(m::Function,c::Function,a::Function,b::Function,
   trial,test)
-  res(t,(u,ut,utt),v) = -1.0 * b(v)
-  jac(t,(u,ut,utt),du,v) = a(du,v)
-  jac_t(t,(u,ut,utt),dut,v) = c(dut,v)
-  jac_tt(t,(u,ut,utt),dutt,v) = m(dutt,v)
+  res(t,u,v) = -1.0 * b(v)
+  jac(t,u,du,v) = a(du,v)
+  jac_t(t,u,dut,v) = c(dut,v)
+  jac_tt(t,u,dutt,v) = m(dutt,v)
   assem_t = SparseMatrixAssembler(trial,test)
   trial_t = ∂t(trial)
   trial_tt = ∂t(trial_t)
@@ -159,10 +159,10 @@ end
 
 function TransientConstantMatrixFEOperator(m::Function,c::Function,a::Function,b::Function,
   trial,test)
-  res(t,(u,ut,utt),v) = m(utt,v) + c(ut,v) + a(u,v) - b(t,v)
-  jac(t,(u,ut,utt),du,v) = a(du,v)
-  jac_t(t,(u,ut,utt),dut,v) = c(dut,v)
-  jac_tt(t,(u,ut,utt),dutt,v) = m(dutt,v)
+  res(t,u,v) = m(∂tt(u),v) + c(∂t(u),v) + a(u,v) - b(t,v)
+  jac(t,u,du,v) = a(du,v)
+  jac_t(t,u,dut,v) = c(dut,v)
+  jac_tt(t,u,dutt,v) = m(dutt,v)
   assem_t = SparseMatrixAssembler(trial,test)
   trial_t = ∂t(trial)
   trial_tt = ∂t(trial_t)
@@ -172,10 +172,10 @@ end
 
 function TransientAffineFEOperator(m::Function,c::Function,a::Function,b::Function,
   trial,test)
-  res(t,(u,ut,utt),v) = m(t,utt,v) + c(t,ut,v) + a(t,u,v) - b(t,v)
-  jac(t,(u,ut,utt),du,v) = a(t,du,v)
-  jac_t(t,(u,ut,utt),dut,v) = c(t,dut,v)
-  jac_tt(t,(u,ut,utt),dutt,v) = m(t,dutt,v)
+  res(t,u,v) = m(t,∂tt(u),v) + c(t,∂t(u),v) + a(t,u,v) - b(t,v)
+  jac(t,u,du,v) = a(t,du,v)
+  jac_t(t,u,dut,v) = c(t,dut,v)
+  jac_tt(t,u,dutt,v) = m(t,dutt,v)
   assem_t = SparseMatrixAssembler(trial,test)
   trial_t = ∂t(trial)
   trial_tt = ∂t(trial_t)
@@ -195,24 +195,23 @@ end
 function TransientFEOperator(res::Function,trial,test;order::Integer=1)
   function jac_0(t,x,dx0,dv)
     function res_0(y)
-      x0 = (y,x[2:end]...)
+      x0 = TransientCellField(y,x.derivatives)
       res(t,x0,dv)
     end
-    jacobian(res_0,x[1])
+    jacobian(res_0,x.cellfield)
   end
   jacs = (jac_0,)
-  for i in 2:order+1
+  for i in 1:order
     function jac_i(t,x,dxi,dv)
       function res_i(y)
-        x = (x[1:i-1]...,y,x[i+1:end]...)
-        res(t,x,dv)
+        derivatives = (x.derivatives[1:i-1]...,y,x.derivatives[i+1:end]...)
+        xi = TransientCellField(x.cellfield,derivatives)
+        res(t,xi,dv)
       end
-      jacobian(res_i,x[i])
+      jacobian(res_i,x.derivatives[i])
     end
     jacs = (jacs...,jac_i)
   end
-  # jac(t,u,ut,du,dv) = jacobian(x->res(t,x,ut,dv),u)
-  # jac_t(t,u,ut,dut,dv) = jacobian(xt->res(t,u,xt,dv),ut)
   TransientFEOperator(res,jacs...,trial,test)
 end
 
@@ -233,10 +232,11 @@ get_order(op::TransientFEOperatorFromWeakForm) = op.order
 function allocate_residual(op::TransientFEOperatorFromWeakForm,uh::FEFunction,cache)
   V = get_test(op)
   v = get_fe_basis(V)
-  xh = ()
+  dxh = ()
   for i in 1:get_order(op)+1
-    xh = (xh...,uh)
+    dxh = (dxh...,uh)
   end
+  xh = TransientCellField(uh,dxh)
   vecdata = collect_cell_vector(V,op.res(0.0,xh,v))
   allocate_vector(op.assem_t,vecdata)
 end
@@ -245,7 +245,7 @@ function residual!(
   b::AbstractVector,
   op::TransientFEOperatorFromWeakForm,
   t::Real,
-  xh::Tuple{Vararg{FEFunction}},
+  xh::TransientCellField,#Tuple{Vararg{FEFunction}},
   cache)
   V = get_test(op)
   v = get_fe_basis(V)
@@ -255,10 +255,11 @@ function residual!(
 end
 
 function allocate_jacobian(op::TransientFEOperatorFromWeakForm,uh::FEFunction,cache)
-  xh = ()
+  dxh = ()
   for i in 1:get_order(op)+1
-    xh = (xh...,uh)
+    dxh = (dxh...,uh)
   end
+  xh = TransientCellField(uh,dxh)
   _matdata = ()
   for i in 1:get_order(op)+1
     _matdata = (_matdata...,matdata_jacobian(op,0.0,xh,i,0.0))
@@ -271,7 +272,7 @@ function jacobian!(
   A::AbstractMatrix,
   op::TransientFEOperatorFromWeakForm,
   t::Real,
-  xh::Tuple{Vararg{FEFunction}},
+  xh::TransientCellField,#Tuple{Vararg{FEFunction}},
   i::Integer,
   γᵢ::Real,
   cache)
@@ -284,7 +285,7 @@ function jacobians!(
   A::AbstractMatrix,
   op::TransientFEOperatorFromWeakForm,
   t::Real,
-  xh::Tuple{Vararg{FEFunction}},
+  xh::TransientCellField,#Tuple{Vararg{FEFunction}},
   γ::Tuple{Vararg{Real}},
   cache)
   _matdata = ()
@@ -319,7 +320,7 @@ end
 function matdata_jacobian(
   op::TransientFEOperatorFromWeakForm,
   t::Real,
-  xh::Tuple{Vararg{FEFunction}},
+  xh::TransientCellField,#Tuple{Vararg{FEFunction}},
   i::Integer,
   γᵢ::Real)
   Uh = evaluate(get_trial(op),nothing)
@@ -342,15 +343,16 @@ function test_transient_fe_operator(op::TransientFEOperator,uh)
   @test isa(U0,FESpace)
   r = allocate_residual(op,uh,cache)
   @test isa(r,AbstractVector)
-  residual!(r,op,0.0,(uh,uh),cache)
+  xh = TransientCellField(uh,(uh,))
+  residual!(r,op,0.0,xh,cache)
   @test isa(r,AbstractVector)
   J = allocate_jacobian(op,uh,cache)
   @test isa(J,AbstractMatrix)
-  jacobian!(J,op,0.0,(uh,uh),1,1.0,cache)
+  jacobian!(J,op,0.0,xh,1,1.0,cache)
   @test isa(J,AbstractMatrix)
-  jacobian!(J,op,0.0,(uh,uh),2,1.0,cache)
+  jacobian!(J,op,0.0,xh,2,1.0,cache)
   @test isa(J,AbstractMatrix)
-  jacobians!(J,op,0.0,(uh,uh),(1.0,1.0),cache)
+  jacobians!(J,op,0.0,xh,(1.0,1.0),cache)
   @test isa(J,AbstractMatrix)
   cache = update_cache!(cache,op,0.0)
   true

src/TransientFETools/TransientFETools.jl

@@ -85,10 +85,22 @@ import GridapODEs.ODETools: GenericODESolution
 import Base: iterate
 export test_transient_fe_solution
 
-export Transient2ndOrderFEOperator
+export TransientCellField
+using Gridap.CellData: CellField
+using Gridap.MultiField: MultiFieldCellField
+using Gridap.FESpaces: FEBasis
+import Gridap.CellData: get_data
+import Gridap.CellData: get_triangulation
+import Gridap.CellData: DomainStyle
+import Gridap.CellData: gradient
+import Gridap.CellData: ∇∇
+import Gridap.CellData: change_domain
+import Gridap.FESpaces: BasisStyle
 
 include("TransientFESpaces.jl")
 
+include("TransientCellField.jl")
+
 include("TransientFEOperators.jl")
 
 include("ODEOperatorInterfaces.jl")

test/DiffEqsWrappersTests/DiffEqsTests.jl

@@ -40,9 +40,9 @@ function fe_problem(u, n)
   b(v, t) = ∫( v * f(t) )dΩ
   m(u, v) = ∫( v * u )dΩ
 
-  res(t, (u, ut), v) = a(u, v) + m(ut, v) - b(v, t)
-  jac(t, (u, ut), du, v) = a(du, v)
-  jac_t(t, (u, ut), dut, v) = m(dut, v)
+  res(t, u, v) = a(u, v) + m(∂t(u), v) - b(v, t)
+  jac(t, u, du, v) = a(du, v)
+  jac_t(t, u, dut, v) = m(dut, v)
 
   op = TransientFEOperator(res, jac, jac_t, U, V0)
 

test/TransientFEsTests/BoundaryHeatEquationTests.jl

@@ -54,9 +54,9 @@ b(v,t) = ∫(v*f(t))dΩ
 m(ut,v) = ∫(ut*v)dΩ
 b_Γ(v,t) = ∫(v*(∇(u(t))⋅nb))dΓ
 
-res(t,(u,ut),v) = a(u,v) + m(ut,v) - b(v,t) - b_Γ(v,t)
-jac(t,(u,ut),du,v) = a(du,v)
-jac_t(t,(u,ut),dut,v) = m(dut,v)
+res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t) - b_Γ(v,t)
+jac(t,u,du,v) = a(du,v)
+jac_t(t,u,dut,v) = m(dut,v)
 
 op = TransientFEOperator(res,jac,jac_t,U,V0)
 

test/TransientFEsTests/DGHeatEquationTests.jl

@@ -54,9 +54,9 @@ b_Γ(v,t) = ∫( (γ/h)*v*u(t) - (∇(v)⋅nb)*u(t) )dΓ
 a_Λ(u,v) = ∫( (γ/h)*jump(v*ns)⊙jump(u*ns) - jump(v*ns)⊙mean(∇(u)) - mean(∇(v))⊙jump(u*ns) )dΛ
 
 
-res(t,(u,ut),v) = a(u,v) + m(ut,v) + a_Γ(u,v) + a_Λ(u,v) - b(v,t) - b_Γ(v,t)
-jac(t,(u,ut),du,v) = a(du,v) + a_Γ(du,v) + a_Λ(du,v)
-jac_t(t,(u,ut),dut,v) = m(dut,v)
+res(t,u,v) = a(u,v) + m(∂t(u),v) + a_Γ(u,v) + a_Λ(u,v) - b(v,t) - b_Γ(v,t)
+jac(t,u,du,v) = a(du,v) + a_Γ(du,v) + a_Λ(du,v)
+jac_t(t,u,dut,v) = m(dut,v)
 
 op = TransientFEOperator(res,jac,jac_t,U,V0)