Transient distributed cell field

.gitignore

@@ -9,3 +9,4 @@
 /tmp/
 *.vtu
 *.pvtu
+*.pvd

Project.toml

@@ -6,6 +6,7 @@ version = "0.2.5"
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
+GridapODEs = "55e38337-5b6e-4c7c-9cfc-e00dd49102e6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"

src/GridapDistributed.jl

@@ -41,4 +41,6 @@ include("DivConformingFESpaces.jl")
 
 include("MultiField.jl")
 
+include("GridapODEs/GridapODEs.jl")
+
 end # module

src/GridapODEs/GridapODEs.jl

@@ -0,0 +1,35 @@
+module GridapODEs
+
+using GridapDistributed: DistributedCellDatum, DistributedCellField, DistributedMeasure
+using GridapDistributed: DistributedSingleFieldFEFunction, DistributedMultiFieldFEFunction
+using GridapDistributed: DistributedSingleFieldFESpace, DistributedMultiFieldFESpace
+using GridapDistributed: local_views
+using Gridap.Fields
+using Gridap.Geometry
+using Gridap.Arrays
+using Gridap.CellData
+using Gridap.Algebra
+using Gridap.FESpaces
+using Gridap.MultiField: num_fields
+using Gridap.Helpers: first_and_tail
+using GridapODEs.TransientFETools
+using GridapODEs.ODETools
+using PartitionedArrays
+
+import GridapODEs.ODETools: ∂t, ∂tt
+import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
+import LinearAlgebra: det, tr, cross, dot, ⋅
+import Base: inv, abs, abs2, *, +, -, /, adjoint, transpose, real, imag, conj
+import Gridap.Fields: grad2curl
+
+export TransientDistributedCellField
+export TransientSingleFieldDistributedCellField
+export TransientMultiFieldDistributedCellField
+
+include("TransientDistributedCellField.jl")
+
+include("TransientMultiFieldDistributedCellField.jl")
+
+include("TransientFESpaces.jl")
+
+end

src/GridapODEs/TransientDistributedCellField.jl

@@ -0,0 +1,99 @@
+# Transient Distributed CellField
+abstract type TransientDistributedCellField <: DistributedCellDatum  end
+
+# Transient SingleField
+struct TransientSingleFieldDistributedCellField{A} <: TransientDistributedCellField
+  cellfield::A
+  derivatives::Tuple
+end
+
+# Constructors
+function TransientFETools.TransientCellField(single_field::DistributedSingleFieldFEFunction,derivatives::Tuple)
+  TransientSingleFieldDistributedCellField(single_field,derivatives)
+end
+
+function TransientFETools.TransientCellField(single_field::DistributedCellField,derivatives::Tuple)
+TransientSingleFieldDistributedCellField(single_field,derivatives)
+end
+
+# Time derivative
+function ∂t(f::TransientDistributedCellField)
+  cellfield, derivatives = first_and_tail(f.derivatives)
+  TransientCellField(cellfield,derivatives)
+end
+
+∂tt(f::TransientDistributedCellField) = ∂t(∂t(f))
+
+# Integration related
+function Fields.integrate(f::TransientDistributedCellField,b::DistributedMeasure)
+  integrate(f.cellfield,b)
+end
+
+# Differential Operations
+Fields.gradient(f::TransientDistributedCellField) = gradient(f.cellfield)
+Fields.∇∇(f::TransientDistributedCellField) = ∇∇(f.cellfield)
+
+# Unary ops
+for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj)
+  @eval begin
+    ($op)(a::TransientDistributedCellField) = ($op)(a.cellfield)
+  end
+end
+
+# Binary ops
+for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/)
+  @eval begin
+    ($op)(a::TransientDistributedCellField,b::TransientDistributedCellField) = ($op)(a.cellfield,b.cellfield)
+    ($op)(a::TransientDistributedCellField,b::DistributedCellField) = ($op)(a.cellfield,b)
+    ($op)(a::DistributedCellField,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
+    ($op)(a::TransientDistributedCellField,b::Number) = ($op)(a.cellfield,b)
+    ($op)(a::Number,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
+    ($op)(a::TransientDistributedCellField,b::Function) = ($op)(a.cellfield,b)
+    ($op)(a::Function,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
+  end
+end
+
+Base.broadcasted(f,a::TransientDistributedCellField,b::TransientDistributedCellField) = broadcasted(f,a.cellfield,b.cellfield)
+Base.broadcasted(f,a::TransientDistributedCellField,b::DistributedCellField) = broadcasted(f,a.cellfield,b)
+Base.broadcasted(f,a::DistributedCellField,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
+Base.broadcasted(f,a::Number,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
+Base.broadcasted(f,a::TransientDistributedCellField,b::Number) = broadcasted(f,a.cellfield,b)
+Base.broadcasted(f,a::Function,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
+Base.broadcasted(f,a::TransientDistributedCellField,b::Function) = broadcasted(f,a.cellfield,b)
+Base.broadcasted(a::typeof(*),b::typeof(∇),f::TransientDistributedCellField) = broadcasted(a,b,f.cellfield)
+Base.broadcasted(a::typeof(*),s::Fields.ShiftedNabla,f::TransientDistributedCellField) = broadcasted(a,s,f.cellfield)
+
+dot(::typeof(∇),f::TransientDistributedCellField) = dot(∇,f.cellfield)
+outer(::typeof(∇),f::TransientDistributedCellField) = outer(∇,f.cellfield)
+outer(f::TransientDistributedCellField,::typeof(∇)) = outer(f.cellfield,∇)
+cross(::typeof(∇),f::TransientDistributedCellField) = cross(∇,f.cellfield)
+
+# Skeleton related
+function Base.getproperty(f::TransientDistributedCellField, sym::Symbol)
+  if sym in (:⁺,:plus,:⁻, :minus)
+    derivatives = ()
+    cellfield = DistributedCellField(f.cellfield,sym)
+    for iderivative in f.derivatives
+      derivatives = (derivatives...,DistributedCellField(iderivative))
+    end
+    return TransientSingleFieldCellField(cellfield,derivatives)
+  else
+    return getfield(f, sym)
+  end
+end
+
+Base.propertynames(x::TransientDistributedCellField, private::Bool=false) = propertynames(x.cellfield, private)
+
+for op in (:outer,:*,:dot)
+  @eval begin
+    ($op)(a::TransientDistributedCellField,b::SkeletonPair{<:DistributedCellField}) = ($op)(a.cellfield,b)
+    ($op)(a::SkeletonPair{<:DistributedCellField},b::TransientDistributedCellField) = ($op)(a,b.cellfield)
+  end
+end
+
+Arrays.evaluate!(cache,k::Operation,a::TransientDistributedCellField,b::SkeletonPair{<:DistributedCellField}) = evaluate!(cache,k,a.cellfield,b)
+
+Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:DistributedCellField},b::TransientDistributedCellField) = evaluate!(cache,k,a,b.cellfield)
+
+CellData.jump(a::TransientDistributedCellField) = jump(a.cellfield)
+CellData.mean(a::TransientDistributedCellField) = mean(a.cellfield)

src/GridapODEs/TransientFESpaces.jl

@@ -0,0 +1,66 @@
+# Functions for transient FE spaces
+
+Fields.evaluate(U::DistributedSingleFieldFESpace,t::Nothing) = U
+
+(U::DistributedSingleFieldFESpace)(t) = U
+
+∂t(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U)
+∂t(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂t.(U.field_fe_space))
+∂tt(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U)
+∂tt(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.field_fe_spaces))
+
+function TransientMultiFieldFESpace(spaces::Vector{<:DistributedSingleFieldFESpace})
+  MultiFieldFESpace(spaces)
+end
+
+# Functions for transient FE Functions
+
+function ODETools.allocate_jacobian(
+  op::TransientFETools.TransientFEOperatorFromWeakForm,
+  duh::DistributedCellField,
+  cache)
+  _matdata_jacobians = TransientFETools.fill_initial_jacobians(op,duh)
+  matdata = _vcat_distributed_matdata(_matdata_jacobians)
+  allocate_matrix(op.assem_t,matdata)
+end
+
+function ODETools.allocate_jacobian(
+  op::TransientFETools.TransientFEOperatorFromWeakForm,
+  duh::DistributedMultiFieldFEFunction,
+  cache)
+  _matdata_jacobians = TransientFETools.fill_initial_jacobians(op,duh)
+  matdata = _vcat_distributed_matdata(_matdata_jacobians)
+  allocate_matrix(op.assem_t,matdata)
+end
+
+function ODETools.jacobians!(
+  A::AbstractMatrix,
+  op::TransientFETools.TransientFEOperatorFromWeakForm,
+  t::Real,
+  xh::TransientDistributedCellField,
+  γ::Tuple{Vararg{Real}},
+  cache)
+  _matdata_jacobians = TransientFETools.fill_jacobians(op,t,xh,γ)
+  matdata = _vcat_distributed_matdata(_matdata_jacobians)
+  assemble_matrix_add!(A,op.assem_t, matdata)
+  A
+end
+
+function _vcat_distributed_matdata(_matdata)
+  term_to_cellmat = map_parts(a->a[1],local_views(_matdata[1]))
+  term_to_cellidsrows = map_parts(a->a[2],local_views(_matdata[1]))
+  term_to_cellidscols = map_parts(a->a[3],local_views(_matdata[1]))
+  for j in 2:length(_matdata)
+    term_to_cellmat_j = map_parts(a->a[1],local_views(_matdata[j]))
+    term_to_cellidsrows_j = map_parts(a->a[2],local_views(_matdata[j]))
+    term_to_cellidscols_j = map_parts(a->a[3],local_views(_matdata[j]))
+    term_to_cellmat = map_parts((a,b)->vcat(a,b),local_views(term_to_cellmat),local_views(term_to_cellmat_j))
+    term_to_cellidsrows = map_parts((a,b)->vcat(a,b),local_views(term_to_cellidsrows),local_views(term_to_cellidsrows_j))
+    term_to_cellidscols = map_parts((a,b)->vcat(a,b),local_views(term_to_cellidscols),local_views(term_to_cellidscols_j))
+  end
+  map_parts( (a,b,c) -> (a,b,c),
+    local_views(term_to_cellmat),
+    local_views(term_to_cellidsrows),
+    local_views(term_to_cellidscols)
+  )
+end

src/GridapODEs/TransientMultiFieldDistributedCellField.jl

@@ -0,0 +1,61 @@
+# Transient MultiField
+struct TransientMultiFieldDistributedCellField{A} <: TransientDistributedCellField
+  cellfield::A
+  derivatives::Tuple
+  transient_single_fields::Vector{<:TransientDistributedCellField} # used to iterate
+end
+
+# Constructors
+function TransientFETools.TransientCellField(multi_field::DistributedMultiFieldFEFunction,derivatives::Tuple)
+  transient_single_fields = _to_transient_single_distributed_fields(multi_field,derivatives)
+  TransientMultiFieldDistributedCellField(multi_field,derivatives,transient_single_fields)
+end
+
+# Get single index
+function Base.getindex(f::TransientMultiFieldDistributedCellField,ifield::Integer)
+  single_field = f.cellfield[ifield]
+  single_derivatives = ()
+  for ifield_derivatives in f.derivatives
+    single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield))
+  end
+  TransientSingleFieldDistributedCellField(single_field,single_derivatives)
+end
+
+# Get multiple indices
+function Base.getindex(f::TransientMultiFieldDistributedCellField,indices::Vector{<:Int})
+  cellfield = DistributedMultiFieldCellField(f.cellfield[indices],DomainStyle(f.cellfield))
+  derivatives = ()
+  for derivative in f.derivatives
+    derivatives = (derivatives...,DistributedMultiFieldCellField(derivative[indices],DomainStyle(derivative)))
+  end
+  transient_single_fields = _to_transient_single_distributed_fields(cellfield,derivatives)
+  TransientMultiFieldDistributedCellField(cellfield,derivatives,transient_single_fields)
+end
+
+function _to_transient_single_distributed_fields(multi_field,derivatives)
+  transient_single_fields = TransientDistributedCellField[]
+  for ifield in 1:num_fields(multi_field)
+    single_field = multi_field[ifield]
+    single_derivatives = ()
+    for ifield_derivatives in derivatives
+      single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield))
+    end
+    transient_single_field = TransientSingleFieldDistributedCellField(single_field,single_derivatives)
+    push!(transient_single_fields,transient_single_field)
+  end
+  transient_single_fields
+end
+
+# Iterate functions
+Base.iterate(f::TransientMultiFieldDistributedCellField)  = iterate(f.transient_single_fields)
+Base.iterate(f::TransientMultiFieldDistributedCellField,state)  = iterate(f.transient_single_fields,state)
+
+# Time derivative
+function ∂t(f::TransientMultiFieldDistributedCellField)
+  cellfield, derivatives = first_and_tail(f.derivatives)
+  transient_single_field_derivatives = TransientDistributedCellField[]
+  for transient_single_field in f.transient_single_fields
+    push!(transient_single_field_derivatives,∂t(transient_single_field))
+  end
+  TransientMultiFieldDistributedCellField(cellfield,derivatives,transient_single_field_derivatives)
+end

test/GridapODEs/HeatEquationTests.jl

@@ -0,0 +1,77 @@
+module HeatEquationTests
+
+using Gridap
+using GridapODEs
+using GridapDistributed
+using PartitionedArrays
+using Test
+
+function main(parts)
+  θ = 0.2
+
+  u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
+  u(t::Real) = x -> u(x,t)
+  f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)
+
+  domain = (0,1,0,1)
+  partition = (4,4)
+  model = CartesianDiscreteModel(parts,domain,partition)
+
+  order = 2
+
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  V0 = FESpace(
+    model,
+    reffe,
+    conformity=:H1,
+    dirichlet_tags="boundary"
+  )
+  U = TransientTrialFESpace(V0,u)
+
+  Ω = Triangulation(model)
+  degree = 2*order
+  dΩ = Measure(Ω,degree)
+
+  #
+  m(u,v) = ∫(u*v)dΩ
+  a(u,v) = ∫(∇(v)⋅∇(u))dΩ
+  b(t,v) = ∫(v*f(t))dΩ
+
+  res(t,u,v) = a(u,v) + m(∂t(u),v) - b(t,v)
+  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)
+  op_constant = TransientConstantMatrixFEOperator(m,a,b,U,V0)
+
+  t0 = 0.0
+  tF = 1.0
+  dt = 0.1
+
+  U0 = U(0.0)
+  uh0 = interpolate_everywhere(u(0.0),U0)
+
+  ls = LUSolver()
+  ode_solver = ThetaMethod(ls,dt,θ)
+
+  sol_t = solve(ode_solver,op,uh0,t0,tF)
+  sol_t_const = solve(ode_solver,op_constant,uh0,t0,tF)
+
+  l2(w) = w*w
+
+  tol = 1.0e-6
+
+  for (uh_tn, tn) in sol_t
+    e = u(tn) - uh_tn
+    el2 = sqrt(sum( ∫(l2(e))dΩ ))
+    @test el2 < tol
+  end
+
+  for (uh_tn, tn) in sol_t_const
+    e = u(tn) - uh_tn
+    el2 = sqrt(sum( ∫(l2(e))dΩ ))
+    @test el2 < tol
+  end
+end
+
+end #module